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Chapter  1 


INTRODUCTION 


Normal  mode  theory  determines  exact  solutions  for  ducted  sound  propagation 
problems.  The  theory  gives  a  characteristic  equation  for  a  given  duct  system. 

Unfortunately,  the  characteristic  equations  which  determine  the  eigenvalues  are  usually 
transcendental  functions  that  require  numerical  analysis  and  the  aid  of  a  ccm^uter  to  solve. 
Using  approximations  of  a  characteristic  equation  can  make  finding  the  eigenvalues  easier 
and  quicker. 

An  important  duct  problem  is  a  three-layer  model  of  the  ocean.  The  characteristic 
equation  determining  its  eigenvalues  is  complex  and  transcendental.  It  can  be  solved 
numerically  using  Newton's  Method  for  nonlinear  systems,  but  that  is  a  time-consuming 
computer  process.  This  study  describes  a  factoring  technique  for  breaking  the 
characteristic  equation  into  simpler  equations,  and  compares  the  eigenvalues  produced  by 
the  characteristic  equation  and  the  simpler  equations. 

Two  types  of  three  layer  models  will  be  considered:  a  liquid/liquid/liquid  (LLL) 
model  and  a  liquid/elastic/elastic  (I  .EE)  model.  The  LLL  nxxlel  is  the  simpler  since  only 
pressure  waves  can  propagate  in  liquids.  The  LEE  treats  the  water  as  a  liquid,  but  treats  the 
sediment  and  the  bedrock  as  elastic  media  which  support  both  pressure  and  shear  waves. 
The  use  for  two  models  arises  from  many  types  of  sediments  and  ocean  bottoms. 
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1.1.  Definition  of  the  Three-Layer  Model 

The  model  representing  the  ocean  is  a  three-layer  waveguide  composed  of  water, 
sediment,  and  bedrock.  The  water  is  designated  as  layer  1  and  has  constants  density  pi, 
sound  speed  ci,  and  wavenumber  ki.  Likewise,  the  sediment  is  layer  2  with  p2,  C2,  and 
k2,  and  the  bedrock,  usually  basalt  or  granite,  is  layer  3  with  pi,  C3,  and  k^.  The  water  has 
depth  h  and  the  sediment  has  depth  t.  The  bedrock  is  a  halfspace.  All  boundary  surfaces 
are  smooth  and  horizontal,  and  the  water-air  surface  at  z=0  is  a  pressure  release  surface. 


air 


water 

Px 

sediment 

Pi 

^2 

9 

9 

9 

9 

bedrock 

P3 

^3 

r  =  0 


z  =  /i 


z  =  h+t 


z 


Figure  1.1.  The  Three-Layer  Ocean  Model 
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1.2.  Development  of  the  Characteristic  Equation 

The  development  of  the  characteristic  equation  for  a  general  layered  waveguide 
begins  with  the  wave  equation's  associated  Helmholtz  equation  for  the  pressure  of  plane 
waves.  Pressure  P  is  a  function  of  both  range  r  and  depth  z.  The  wavenumber  k  and 
sound  speed  c  are  generally  functions  of  z,  and  k(z)=oilc(z).  Since  the  interest  is  in  the 
water  layer,  the  equation  is  developed  for  this  region. 

vV(r.z)  +A:i(z)2p(r,z)=0  .  (1.1) 

The  separation  of  variables  technique  can  be  used  to  get  the  differential  equation  for  Z(z) 
that  describes  the  z-dependence  only  as  in  (1.2).  The  separation  constant  is  denoted 
and  is  the  horizontal  component  of  the  wavenumber  k.  The  kn^  are  the  eigenvalues  that 
determine  the  normal  modes  for  the  waveguide.  For  example,  in  the  water  layer, 

+  iki{zf-khZi{z)  =  0  .  (1.2) 

dz^ 

The  vertical  component  of  ki  will  be  denoted  as  JQ.  The  relationship  for  k  and  its 
components  is  ki  ^=Ki  ^+kn^.  For  the  models  developed  here,  k  and  ^rwill  not  be 
functions  of  z,  since  the  sound  speeds  are  assumed  to  be  constants. 

^^i^+KfZi(z)  =  0  .  (1.3) 

dz^ 

The  solution  to  (1.3)  is  given  by  (1.4),  where  represents  an  upward 
traveling  wave,  and  represents  a  downward  traveling  wave. 


Zi(z)=i4e‘*»*  +  Be-*'f«*  . 


(1.4) 
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If  the  water  region  is  divided  into  a  lower  and  upper  region  at  the  depth  of  a  source,  then 
two  equations  can  be  written  from  (1.4):  one  for  a  wave  that  travels  upward  from  the 
source  and  strikes  the  surface,  and  one  that  travels  downward  from  the  source  and  strikes 
the  water-sediment  interface.  The  uptraveling  wave  becomes  a  downward  traveling  wave 
after  reflecting  off  the  surface  and  the  down  traveling  wave  becomes  an  uptraveling  wave 
after  hitting  the  interface.  The  two  equations  describing  this  are 

=  ,  (1.5) 


and 


»>n(z-A  )  +  g-»ri(z-A  ) 


(1.6) 


where  SRi  and  9ti3  are  complex  reflection  coefflcients.  To  relate  these  two  equations  to 
one  another,  their  Wronskian  is  set  equal  to  zero  in  (1,7). 


W(/Ci)  = 


Z^' 


z^' 


=  Z^Z^'-Z^'Z^  =  0 


(1.7) 


After  completing  the  necessary  math,  the  characteristic  equation  is 

W(K'i)  =  =  0  .  (1.8) 


1.2.1.  Making  kn  Complex 

Now  kn ,  the  horizontal  component  of  wavenumber  k,  is  made  complex  with  small 
imaginaiy  part  Allowing  the  wavenumber  to  be  complex  introduces  a  loss  into  the 
characteristic  equation. 


kn  —  kn  ^  tCn  . 


(1.9) 
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The  expression  for  the  vertical  component  of  the  wavenumber  is  then 

.  (1.10) 

If  it  is  assumed  that  £V,  is  very  small  so  its  squared  term  may  be  ignored,  after  rearranging 
and  using  the  binomial  expansion,  the  expression  for  x*!  becomes 

Kx={k\-k})^^-'-^  .  (1.11) 

The  first  term  on  the  right  side  of  (1.1 1)  represents  the  real  part  of  Ki  and  will  be  denoted 
as  xf .  The  Ki  in  the  denominator  of  the  second  term  on  the  right  side  of  ( 1 . 1 1)  can  be 
approximated  as  xf  so  that  rewriting  (1.11)  gives 

V,  =  xf-^^  .  (1.12) 

Substituting  (1.12)  into  the  characteristic  equation  gives 

=  0  .  (1.13) 

The  exponential  with  the  £n  in  it  is  the  loss  ot  attenuation  introduced  by  the  complex  kn. 
The  modes  with  the  smallest  £^'s  will  propagate  best. 

1.2.2.  The  Magnitude  and  Phase  of 

The  reflection  coefficient  SRi  is  complex  and  may  be  represented  in  polar  form. 

=  .  (1.14) 

The  air-water  interface  is  considered  to  be  a  near  perfect  pressure  release  boundary  for  all 


calculations.  Therefore  the  value  of  SHi  is  set  equal  to  0.99  and  the  phase  is  set  equal  to 
pi.  The  characteristic  equation  is  now 

1 -|9li|9?i3e-‘V+^i)e-2*-WK^  =  0  .  (1.15) 

1.3.  Solving  the  Exact  Characteristic  Equation 

Before  the  characteristic  equation  can  be  solved,  SR13  must  be  known.  9{i3  will  be 
developed  later,  and  will  be  a  complex  function  of  both  k„  and  e„.  The  characteristic 
equation,  once  5Ri3  is  substituted  into  it,  can  be  separated  into  its  real  and  imaginary  parts 
yielding  two  equations  for  the  two  unknowns,  and  £^.  Newton's  Method  for  nonlinear 
systems  (Burden  and  Faires  1985,  p.  496}  can  be  used  to  solve  for  kn  and  £^. 

Recall  that  for  a  linear  system,  Newton's  Method  is  an  iterative  process  in  which 
Xi's  are  sought  that  satisfy 

/(x.)  =  0,  (1.16) 

,  (1.17) 

f'Oci) 

(Kunz  1957,  pp.  11),  where  Xi  is  the  initial  guess  and  Xi^i  becomes  the  new  guess  after 
each  iteration.  For  the  two  variable  case,  two  equations  are  needed. 


flikntCn)  ~  0 

(1.18) 

(1.19) 

The  iterative  equatitxi  in  matrix  notation  is 


-1 

^ni*\ 

kn 

dkn  den 

.  ^+1  . 

.£n. 

dft  9/2 

1/2  J 

.  dkn  d£n  . 

(1.20) 


Defining  x  as  in  (1.21),  (1.20)  can  be  rewritten  as  (1.22). 


X  = 


n 

I  Bn 


(1.21) 


G(x)  =  x-J-»(x)F(x)  .  (1.22) 

The  matrix  J(x)  is  the  Jacobian  matrix,  and  J*l(x)  is  the  inverse  of  the  Jacobian.  This 
method  usually  gives  quadratic  convergence  provided  a  sufficiently  accurate  starting  value 
is  given  and  J'^x)  exists  (Burden  and  Faires  1985,  pp.  498). 

1.4.  The  Factoring  Technique 

In  Chapters  2  and  3, 5ii3  will  be  shown  to  be  the  sum  of  two  or  more  terms.  For 
simplicity,  say  9li3  is  equal  to  x  plus  y,  and  and  are  equal  to  unity.  Then  the 

general  form  of  the  characteristic  equation  is 

l-x~y  =  0  .  (1.23) 

If  we  assume  that  x  and  y  are  small,  and  that  the  product  of  the  two  is  also  small,  then  the 
product  xy  can  be  added  to  the  left  side  of  (1.23).  The  right  side  is  kept  at  zero, 
introducing  product  error  a^.  The  left  side  can  be  factored  into  (l-x)(l-y).  Adding  xy  to 
the  left  side  and  then  factoring  in  this  fashion  is  the  factcxing  technique  used  on  the 
charactCTistic  equation.  Each  of  the  factors  (1-x)  and  (1-y)  is  then  separated  into  its  real 
and  imaginaiy  parts.  The  equation  frcm  taking  the  imaginary  part  is  always  a  function  of 
kn  only,  as  will  be  shown  later,  and  can  be  solved  using  Newton's  linear  method.  The 


8 


values  of  kn  can  then  be  substituted  into  the  real  part  equation,  which  is  a  function  of  the 
two  variables,  to  find  the  corresponding  values  of 


9 


Chapter  2 


THE  LIQUID/LIQUID/LIQUID  MODEL 


In  the  liquid/liquid/liquid  (LLL)  model,  all  three  layers  in  the  waveguide  are  liquids 
capable  of  supporting  the  propagation  of  compressional  waves.  The  reflection  and 
transmission  coefficients  needed  in  SR13  in  the  characteristic  equation  must  therefore  be  for 
waves  striking  liquid/liquid  boundaries.  Shear  waves  cannot  propagate  in  liquids  so  they 
are  not  considered  until  Qianter  3. 

2.1.  The  LLL  Characteristic  Equation 

In  order  to  solve  the  characteristic  equation  developed  in  Chapter  1  (1.8),  9li3  must 
be  obtained.  9Ii3  is  a  combined  reflectitm  coefficient  describing  the  total  returned  sound 
from  the  sediment  and  basement  The  derivation  follows  fix>m  Clay  and  Medwin  (1977, 
pp.  96). 

2.1.1.  9^13  for  the  LLL  Model 

The  reflection  coefficient  and  tran^ssion  coefficient  for  two  liquid  layers  are 

Z2-Z1 
Z2  +  Z1 


91i2  = 


(2.1) 
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3i2  = 


2Z2 

Z2  +  Zi 


(2.2) 


where  the  impedances,  Zi  and  Z2,  are  giver  in  Table  2.1.  The  energy  equation  that  relates 
these  coefficients  is  also  listed  in  Table  2.1.  The  phase  that  arises  from  reflection  off  the 
water-sediment  interface  ^2  is  assumed  to  be  zero.  Similarly,  the  equations  for  9I23  and 
323  and  their  energy  relationship  are  in  Table  2.1;  023  is  also  assumed  zero.  As  shown  in 
Figure  2,1,  the  angles  are  measured  from  the  horizontal,  and  the  sediment  compressional 
wave  path  distance  in  between  reflections  or  transmissions  is  r^. 

From  Figure  2.1,  the  total  up  traveling  signal  is  the  sum  of  an  infinite  number  of 
reflections  and  transmissions.  Each  up  and  each  down  path  in  the  sediment  has  a  phase 
delay  of  2k2smd2t  or  2>C2/.  Assuming  the  incident  wave  has  unit  amplitude,  the  total 
reflection  9ii3  is 

91i3  =  91i2  +  S129I23321  ^  +  3i29^23^2i321  •  (2.3) 

After  factoring  out  3 12^23321  e  from  all  but  the  first  term,  the  remaining  terms  in 
(2.3)  have  the  form  of  a  geometric  series  (Spiegel  1968,  pp.  107). 

£  jc-^  =  (l-x)-^  .  (2.4) 

)=o 

The  complex  coefticients  can  be  replaced  by  their  magnitudes  since  all  of  their 
phases  are  approximately  by  zero.  Rewriting  (2.3)  by  summing  the  series  gives  for  91 13 

9i.3=9t.2^  3.2  321 9123  (2  5) 

1  —  9I21  9^23  ^ 


The  reflection  coefficient  is  an  oscillating  functitm  that  depends  on  iKzt.  It  also  depends  on 
the  frequency  and  angle  of  incidence  for  a  given  layer  (Clay  and  Medwin  1977,  pp.  67). 


Table  2.1. 
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For  the  liquid/liquid  boundary  case,  9l2i  =  “  ^12  and  3i232i  =  1-91 12^- 
Substitution  into  (2.5)  yields 


91i3  =  91i2  + 


(I-9I12  )9l23g-^*^ 
1  +  91i2  9I23  g 


(2.6) 


2.1.2.  Incorporation  of  Attenuation  into  9ti3  for  LLL  Model 


Whenever  a  signal  travels  through  the  sediment  it  is  attenuated.  The  values  for  the 
sediment  attenuation  coefficient  Oz,  usually  given  in  dB/(m  kHz),  vary  for  different  types 
of  sediments.  An  attenuation  factor  will  be  necessary  in  the  same  places  as  the  sediment 
phase  factor  in  the  equation  for  S^is. 
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^  _  lQ-i‘*ih0)(fk000)r2 

(2.7) 

a?  _  o?  .  (1  -  9(12^) 

Vll3  =  5ii2  + - - -  . 

1  +  91i2  9123-A  e“2i>f2/ 

(2.8) 

Due  to  the  complex  kr,  xi  is  replaced  by 

=  ^  ■ 

(2.9) 

Making  the  following  simplifications  transfonns  (2.8)  into  (2.12). 

5,  =/l^l-S,2^)3i23  . 

(2.10) 

^2  =A^9l|2  9l23  , 

(2.11) 

<»  CB  5, 

J'u  =  91i2  + - ‘ ; - - —  . 

l  +  S2e-2‘'cfte-2W/Kf 

(2.12) 

2.1.3.  Substitution  of  91 13  into  the  LLL  Characteristic  Equation 

Before  9li3  can  be  substituted  into  the  characteristic  equation,  the  second  term  is 
multiplied  by  the  complex  ccmjugate  of  its  denominator. 

5li3  =  9li2  +  5i£|(«-2‘'^'  +  52£|)/ben  .  (2.13) 

where  £2  and  Den  are 

=  ,  (2.14) 


Den  =  1  +  252^2  cos(2jiC20  +  S2  £2  • 


(2.15) 
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Substituting  (2.13)  into  the  characteristic  equation  gives 

1  -  =  0  .  (2.16) 

Den 


where  E\  and  £12  are 


=  ,  (2.17) 

£12=  £f£^(e-‘(2‘f^*+2Kf»i-<Jh)+52£|e-«(2>c^*+2<*h))  .  (2.18) 


2.2.  Solving  the  Characteristic  Equation 

Equation  (2.16)  is  a  complex  transcendental  equation  whose  solutions  require  the 
Newton’s  method  for  two  variables  described  in  Chapter  1.  The  solutions  found  by  this 
Newton’s  method  will  be  called  the  “exact”  solutions.  To  get  two  equations  for  F(x),  the 
characteristic  equation  is  separated  into  its  real  and  imaginary  parts.  The  real  part  equation 
from  the  characteristic  equation  is 

FI  =  cos(2A:i+0i) 

2  (2.19) 

-  [cos(2^,+2Ar2+<l>i)  +  52£|cos(2^:]+<I>i)]  =  0  , 

Den 

where  K\  and  Kt  are 

ATi  =  /rf  A  ,  (2.20) 


K2-  K^t  , 


(2.21) 
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and  the  imaginary  part  equation  is 

F2  =  SP.12  sin(2A:i+d>i)+  ^^[sin(2A'i+2/:2+<l>i)+  SiEI  sin(2^i+0i)]=  0  .  (2.22) 

Den 

The  partial  derivatives  involved  in  finding  the  Jacobian  are  difficult  to  obtain,  since 
the  xj’s  are  functions  of  kn  and  because  of  the  many  products.  Since  the  Jacobian  is  just  a 
focusing  mechanism  for  determining  the  next  guess,  the  reflection  and  transmission 
coefficients  are  held  constant  with  respect  to  the  it^’s  for  the  derivatives.  The  consequence 
of  this  may  be  more  iterations  necessary  to  focus  on  a  solution. 

Newton's  Method  also  requires  initial  guesses  to  start  the  iterations.  In  the 
program,  the  water  factored  solutions  are  calculated  first,  providing  the  initial  guesses  for 
the  sediment  solutions.  Both  the  water  and  sediment  solutions  are  used  for  initial  guesses 
for  the  exact  solutions.  See  the  Appendix  for  the  program  listing. 

2.3.  Factoring  the  LLL  Characteristic  Equation 

To  factor  (2.16)  as  described  in  section  1.4,  a  cross  term  must  be  added  to  the  left 
hand  side.  Then  the  equation  of  factors  is  (2.23). 

Den 

Each  of  the  factors  can  then  be  independently  set  equal  to  zero.  The  first  factor  is  called  the 
water  factor  since  all  the  parameters  in  it  are  associated  with  the  wave  paths  in  the  water. 

Its  solutions  will  be  called  the  water  solutions.  The  second  factor  will  be  called  the 
sediment  factor  and  its  solutions,  the  sediment  solutiois. 


2.3.1.  The  Water  Factor 


To  solve  for  the  eigenvalues  of  the  water  factor , 

(1 =0  ,  (2.24) 

the  real  and  imaginary  parts  are  set  equal  to  zero. 

1  -  £?  cos(2^i+<l>i)  =  0  .  (2.25) 

sin(2A:i+<t>i)  =  0  .  (2.26) 

The  equation  formed  by  the  imaginary  part  is  a  function  of  kn  only,  which  is  solved  for 
directiy. 

The  value  of  Infcis  chosen  instead  of  n;rin  order  to  make  cos(2/t;r)  always  equal  to  unity 
which  can  be  substituted  into  equation  (2.25)  along  with  kn  to  get 

€„=  In(9li9li2)^  .  (2.28) 


(2n;r-<I>i) 

2h 


(2.27) 
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At  this  point,  another  approximation  is  made  to  make  this  factor  easily  solvable.  The 
variable  S2  is  set  equal  to  zero.  Setting  S2  equal  to  zero  is  equivalent  to  only  keeping  the 
first  bottom  bounce  contribution  in  Figure  2.1..  Equation  (2.29)  becomes 

(1  -  9liSi£f£|c-‘(2^^i+2i:2+<ih)  =  0  .  (2.30) 

Taking  the  real  and  imaginary  pans  of  (2.30)  yields  the  following  two  equations. 

1  -  ^iSiEfElcosi2Ki+2K2+<I>i)  =  0  .  (2.31) 

sin(2^:i+2^2+<l^i)  =  0  .  (2.32) 

The  imaginary  part  equation  is  only  a  function  of  kn  and  can  be  solved  using  Newmn's 
method  for  one  variable  (Kunz  1957,  pp.  11).  Following  the  example  of  the  water  factor, 
the  argument  of  the  sine  function  is  set  equal  u>  2n7c.  Also  xf  and  xf  are  replaced  by  their 
k„  forms  to  get 

Fikn)  =  2h(k^-k^)^K  2r(k|-k2)*/2+  <Pi  -  2n;r  =  0  .  (2.33) 

Newton's  method  requires  the  derivative  of  F  and  initial  guess.  The  initial  guesses  are 
provided  by  the  kn's  found  by  the  water  factor.  See  the  Appendix  for  the  program  listing. 
Once  the  sediment  factor  kn's  are  found,  they  can  be  placed  in  (2.31),  solved  for 

Cn> 

1 _ 


£W=  ln(9liSi) 


(2.34) 
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2.4.  A  Qualitative  Comparison  of  the  Solutions 

The  shapes  of  the  three  types  of  mode  sets,  exact,  water,  and  sediment,  varies. 
Figure  2.2  shows  the  typical  shapes  for  each  set  The  number  of  modes  is  determined  by 
the  frequency/and  the  water  depth  h  by  (2.35)  (Clay  and  Medwin  1977,  pp.  306). 

rt  ^  ^  (2Asin6c)  +  ^  .  (2.35) 

The  modes  are  always  close  together  for  high  jb^’s  and  spread  out  for  low  kn's.  Also,  the 
harder  the  sediment,  the  fewer  the  modes  because  the  minimum  angle  Be  for  which  there  is 

a  mode  increases.  The  discontinuity  in  the  water  solutions  occurs  at  the  kn  which 
corresponds  to  Sc- 

In  Figure  2.2  and  all  the  eigenvalue  plots,  the  higher  up  on  the  y  axis  a  point  is,  the 
better  that  mode  propagates,  i.e.,  the  1^  I  is  smaller  so  the  loss  is  lower  for  that  mode.  The 
range  of  the  Cn  values  is  several  oiders  of  magnitude,  therefore  a  log  scale  is  used.  Note 
also  that  for  the  exact  and  sediment  factor  solutions,  only  eigenvalues  above  the  critical 
angle  are  plotted.  The  critical  angle  shows  the  onset  of  sediment  penetration. 

2.4.1.  Forcing  the  Water  Solutions  to  the  Exact  Solutions 

If  the  sediment  is  very  dense  and  the  sediment  attenuation  is  very  high,  a  signal  will 
tend  to  stay  in  the  water.  Therefore  the  water  factor  solutions  will  propagate  better  than  the 
sediment  factored  solutions.  The  high  density  of  the  sediment  makes  SR  12  nearly  equal  to 
one  and  thus  the  transmission  into  the  sediment  is  small.  Whatever  waves  do  get  into  the 
sediment  will  die  out  quickly  due  to  the  high  attenuation.  Figure  2.3  shows  the  water, 
sediment,  and  exact  solutions  for  these  inputs.  The  matching  of  the  water  solutions  beyond 
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Eigenvalue  Plot  for  a  Medium  Sediment 


0.00  0.05  0.10  0.15  0.20  0.25 

Re(eigenvaiue) 


Figure  2.2,  Typical  T  .1 T .  Eigenvalue  Plot  for  a  Medium  Sediment 


LLL  Forcing  Water  Solutions  to  Exact  Solutions 


0.00  0.02  0.04  0.06  0.08  0.10  0.12  0.14 

Re(eigen  value) 


Figure  2.3.  Forcing  the  LLL  Wato*  Solutions  to  the  Exact  Solutions 
(P2/Pl=7,  C2/Ci=2.33,  012=1) 
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the  critical  angle  with  the  exact  solutions  indicates  that  the  factors  may  be  a  good 
approximation  of  the  exact  solutions. 

2.4.2.  Forcing  the  Sediment  Solutions  to  the  Exact  Solutions 

If  the  sediment  density  and  sound  speed  are  set  approximately  equal  to  the  water 
density  and  sound  speed,  and  a  small  sediment  attenuation  is  chosen,  the  sediment  wave 
path  will  dominate.  Most  of  the  initial  wave  is  transmitted  into  the  sediment  because  SR  12  is 
small  for  these  conditions.  Figure  2.4  shows  the  sediment  and  exact  solutions  to  be  almost 
equal  for  these  inputs,  again  indicating  that  the  factors  might  be  a  good  model  of  the  exact 
solutions. 


LLL  Forcing  Sediment  Solutions  to  Exact  Solutions 


Re(eigenvalue) 


water 

sediment 

exact 


Figure  2.4.  Forcing  the  LIJ.  Sediment  Solutions  to  the  Exact  Solutions 
(P2/Pl=l.  C2/Ci*l,  a2*0.01) 
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Chapter  3 


THE  LIQUID/ELASTIC/ELASTIC  MODEL 


The  LEE  model  treats  both  the  sediment  layer  and  the  basement  halfspace  as  elastic 
materials  capable  of  supporting  shear  waves.  Four  new  parameters  are  thus  introduced,  the 
shear  speeds  of  sound  in  the  sediment  and  basement,  C2s  trss.  and  the  shear  attenuation 
coefficient  for  the  sediment  is  a2S- 

It  may  seem  as  though  this  study  is  incomplete  since  it  does  not  address  a  liquid/ 
liquid/elastic  model.  That  case  was  indeed  studied  but  its  results  proved  to  be  only  slightly 
different  than  the  liquid/liquid/liquid  model. 

3.1.  The  LEE  Characteristic  Equation 

The  characteristic  equation  for  the  LEE  nx)del  is  (1.15).  The  derivation  for  for 
this  case  differs  firom  the  LT.T.  case  because  the  refectitm  and  transmission  coefficients  are 
now  for  liquid/elastic  boundaries.  Also,  new  coefficients  are  needed  to  account  for  the 
shear  waves  generated  from  compiessional  waves  at  the  liquid/elastic  and  elastic/elastic 
interfaces. 

When  a  pressure  wave  strikes  a  boundary  of  an  elastic  medium  from  a  liquid 
medium,  not  only  are  pressure  waves  reflected  firm  and  transmitted  through  the  boundary, 
but  also  a  shear  wave  is  transmitted  (Miklowitz  1978,  pp.  156).  If  the  initial  medium  is 
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also  elastic,  there  will  also  be  a  reflected  shear  wave.  These  two  events  are  illustrated  in 
Figure  3. 1.  Pressure  waves  are  drawn  as  solid  lines,  SV  waves  as  dashed  lines.  The 
notation  S  V  stands  for  a  shear  wave  with  vertical  polarization.  Shear  waves  with 
horizontal  polarization  are  not  coupled  to  the  pressure  waves  and  are  not  considered  here 
(Tolstoy  1973,  pp.  188).  An  incident  shear  wave  generates  pressure  waves  at  the 
boundaries  in  addition  to  reflecting  and  transmitting  shear  waves  (Miklowitz  1978,  pp. 
156).  This  is  illustrated  in  Figure  3.2. 

3.1.1.  Reflection  Coefflcients  for  the  LEE  Model 

For  the  boundary  events,  new  liquid/elastic  and  elastic/elastic  reflection  and 
transmission  coefficients  are  required.  The  coefficients  have  been  obtained  from  three 
sources  and  are  rewritten  here  with  consistent  notation  and  with  incident  angles  measured 
fr-om  the  horizontal.  The  coeffidents  are  listed  in  Table  3.1  (Brekhovsldkh  1980,  pp.  43- 
47)  and  Table  3.2  (Miklowitz  1978,  pp.  160-161).  The  new  coeffidents  are  illustrated  in 
Figure  3.3.  Since  signals  that  penetrate  the  bedrock  do  not  return,  the  323’s  are  not  used  in 
the  calculations  and  therefore  are  not  included  in  Table  3.2.  The  energy  relationships  that 
govern  the  behavior  of  the  liquid/elastic  and  elastic/elastic  coeffidents  are  listed  in  Table  3.3 
(Ergin  1952,  pp.  350)  (Miklowitz  1978, 161-162). 

3.1.2.  SK 13  for  the  LEE  Model 

As  in  the  LLL  model,  the  incident  pressure  wave  follows  the  path  shown  in  Figure 
2.1.  The  total  reflection  coefficient  from  the  wave  is  similar  to  (2.3),  assuming  again  that 
the  phases  ^2  and  ^  are  zero.  The  difference  from  (2.3)  is  that  the  coefficients  are  now 
for  liquid/elastic  or  elastic/elastic  boundaries.  The  waves  which  reflea  <x  transmit 
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Figure  3.3.  Illustration  of  LEE  Reflection  and  Transmission  Coefficients 


pressure  waves  are  designated  pp.  The  contribution  to  for  the  incident  pressure  wave 
which  remains  a  pressure  wave,  the  pp  series,  is  then 


+ 


3i232?9tge-2«^ 
1-9121  3^23 


(3.1) 


In  the  elastic  case,  9121  does  not  equal  -9ti2  and  no  further  simplification  takes  place. 

For  an  incident  pressure  wave  that  transmits  a  shear  wave  which  remains  a  shear 
wave  in  the  sediment  as  shown  in  Figure  3.4,  the  ss  wave,  the  contribution  to  9ti3  is 

_ss_ 


assuming  phases  <l>i2^^  and  are  zero. 
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7  - 

L\  = - 

sin  01 


sin  02 


Z2S  = 


P2C2S 

sin 


=  [Z2  cos2(202s)  +  Z2S  Sin2(2^s)  -  Zi]Id 


SrJJ  =  [  Zi  -  Z2  cos2(2%)  +  Z2S  sin2(202s)  ]/d 


=  -  [  Zi  +  Z2  cos2(2^)  -  Z2S  Sin2(2^)  ]/d 


3'5.[2P,|cos(2%)]/p 


3lf  =  [-2P2^sin(2fe)]/B 


D  =  Z2  cos2(2^)  +  Z2S  sin^(2^)  +  Zj 


Table  3.2. 


Reflection  Coefficients 
for  the 

Elastic/Elastic  Boundary 


9123  =  [(Pi-P3)(<72+<74)  -  (P2+P4)(<7i-<73)]/DD 
9123  =[~(Pl+P3)(^2-^4)  +  (P2-P4)(^  1 +<73)]/£)D 
9123  =  [2(P3<7i-Pi  <73)]/DD 
9lS  =  C2(P4<72-P2  <74)]/DD 
DZ)  =[(P1+P3)(<72+<74)  -  (P2+P4)(<71+<73)] 


P2  =  P2c|s 
P3  =P3c|s 
f)2  =  tan^^s  -  1 
f>3  =  tan^^s  -  1 
P\  =  e  tan^  tan^ 
pz-h  tan^  tan^ 

<7i  =  fimdi  tan02S  tan^s 
<73  =  -g  tan^s 


e  =  2p2  +  P3^3 
/=  2(/t2-P3) 
g  =  P2^2  -  P3^ 
h  =  2/i3  +  |i2^ 

P2  =  « tan^ 

P4  =  -/tan^  tan^  tan^ 

<72  =  A  tan^  tan%s 

<74  =Pi 


Table  3.3. 


Energy  Relationships  for  LEE  Model 


Liquid/Elastic 


Pressure  wave  in  water  against  solid 


,  ft>PP2.  P2tan^  „pp2,  p2tan^s  „ps2 
1  =y\i2  + - Jl2  + - 012 

Pitan^i  Pitan^i 


Pressure  wave  in  solid  against  water 


anBis  I  gPf 


tan^ 


P2tan^ 


Shear  wave  in  solid  against  water 


tan^  P2tan62S 


Elastic^astic 


Pressure  wave  in  sediment  against  basement 

1 = 35f 

tan^  P2tan^  P2tan^ 


Shear  wave  in  sediment  against  basement 

I  =  9i|f\  JBh.  9!^+  P^^L  3g^  Pl^  Sif 
tan^s  P2tan^s  P2tan^s 
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Two  combination  shear/pressure  waves  are  also  included  in  SHn,  the  sp  wave  and 
the  ps  wave.  (There  are  many  other  combinations  possible.  The  pp,  ss,  sp,  and  ps  are  the 
only  ones  studied  here.)  The  sp  wave  travels  as  indicated  in  Figure  3.5.  Assuming  its 
phase  changes  at  the  boundaries  and  arc  zero,  the  sp  contribution  to  9ti3  is 

gjSP  ^  3^132? 

l-'K2l9?23e-2‘'f2^ 

Similarly  from  Figure  3.6  and  setting  and  equal  to  zero,  the  fourth  component 
of  9?i3  is 

gjPS^  3?23i? 

l-9l|?9l23C-2”^s/ 

Summing  the  four  terms  together  gives  for  the  LEE  case 

gj  _  ,  3i232i9C23g~^**^  ^ 

l-9l2f9l23e-2««3^  1- gill  9^23 

(3.5) 

,  3i|32?9t23g-**Q^e-<*2s/ ^  3?23|? 

1  -  9121  9123  1  -  9^1?  9123 


3.1.3.  Incorporation  of  Attenuation  into  9{i3  for  the  LEE  Model 

Just  as  in  the  liquid  layers,  signals  in  elastic  layers  are  attenuated.  The  attenuation 
in  the  sediment  is  dependent  upon  the  length  of  the  path  traveled  and  will  be  included 
everywhere  in  9li3  that  the  phase  changes  for  the  paths  occur.  The  path  distance  traveled  is 
r2s.  The  attenuation  factor  for  the  elastic  sediment  for  a2S  in  dB/m  kHz  is 


As=  1(H“'“/io^Aooo)^»  . 


(3.6) 


Figure  3.4.  Reflection  and  Transmission  Paths  for  the  SS  Series 


Figure  3.5.  Reflection  and  Transmission  Paths  for  the  SP  Series 
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Figure  3.6.  Reflection  and  Transmission  Paths  for  the  PS  Series 


Due  to  the  complex  A:,,  X"2  and  ic'2s  arc  replaced  by  the  following  two  equations. 


K-2=Kj-iWKf.  (3.7) 

xfs  • 

Malting  the  following  definitions  and  multiplying  denominators  by  their  complex 
conjugates  gives  for  9ii3 

5'’'“  =  3i23"9tSM2  ,  (3.9) 

SSS,3«3W9tg^2  _  (3  10) 


SS'’-3;i3H9!S-ti4s  , 


(3.11) 


31 


’  =  3i232i9lMi4As  , 

(3.12) 

(3.13) 

|s=9i|?9lgA|  , 

(3.14) 

£25  =  &</>[&  ^ 

(3.15) 

II 

(3.16) 

II 

(3.17) 

3*13 + 


+ 


5PP(£ig-2*^-5r) 

1  -  E|)cos(2xJr)  +  (5|**£|)^ 

5SS£|s(e-2iA:»-5|S) 


1  -  (Sf £^)cos(2»4f)  +  (5f  £|s)2 
^  5  ^**£2Sg~**^”(£2g~*^^2 

1  -  (5f‘’£f)cos(2j^f)  +  (5j**£|)^ 

^  5  PS£2e-fg»(£2sg-«^”-5|^£^^e«^») 


1  -  (5f £^5)005(245?)  +  (5f £^)^ 


(3.18) 


3.1.4.  Substitution  of  9{i3  into  the  LEE  Characteristic  Equation 


Making  a  few  more  definitions  and  placing  (3.18)  into  (1.15)  gives  the  LEE 
characteristic  equation  (3.21). 


Deni  =  1  -  (5j**£|)cos(2i^f)  +  (S?EI)^  . 


(3.19) 
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Den2  =  1  -  (Sf  £^)cos(2)40  +  (5?^2s)^  •  (3-20) 

(£2£2g_i(2«:,+2#:2-Kli,)  _  SjP£fe-‘(2Jf>+<ft)) 

_  (£2£2^e-i(2is:,+2Ar2s+<ft)  _  sf  £fe-‘(2^:i+<ft))  (3.21) 

^)cnl 

_?^(£2£2£2sg-/(2Ar,+Jir2+A:2s+<lH)  _  5f£2£2£-^e-‘(2A^‘+^*-^»+<«^))  =  0  . 

Two  approaches  were  taken  to  solve  this  complex  function  of  both  kn  and  £^.  First, 
Newton’s  method  for  nonlinear  systems,  and  second,  the  factoring  technique. 

3.2.  Solving  the  LEE  Exact  Equation  with  Newton*s  Nonlinear  Method 

To  solve  the  LEE  characteristic  equation  using  Newton’s  Method  for  nonlinear 
systems  as  described  in  section  1.3,  the  equation  must  be  broken  into  its  real  and  imaginary 
parts  to  give  two  equations  for  the  two  unknowns,  kn  and  £^.  The  Jacobian  for  these  two 
equations  must  also  be  determined.  Since  the  calculations  are  done  by  computer,  a 
discussion  of  some  of  the  computational  pitfalls  is  also  included. 

The  imaginaiy  part  of  (3.21)  gives  for  FI: 

FI  =  0  =  9li9li2£f  sin(2A:i+d>i) 

®  cPP 

+  [£f£|sin(2Ari+2£2+<*>i)  -  sin(2£i+d>i)] 

Dcnl 
9t  cSS 

+  [£f£jssin(2A:i+2A:2s+<*>i)  -  5f  £fsin(2A:i+<l>i)]  (3.22) 

+'^^^[E}E2E2ssiM2Ki+K2+K2s+<I>i)-SrEiE2E2ssini2Ki-K2+K2s+<Pi)] 

Deni 

«,CPS  .  oi.  «  1 

+:2^  [E\E2E2ssini2Ki’^‘K2+K2s+Oi)-S^ElE2E^smi2Ki-^K7rK2S+^i)]  . 
Dcn2 
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and  the  real  part  gives  F2: 

F2  =  0  =  1  -  ^i^nElcos(2Ki  +  <I>i) 

»  cPP 

-  Sr£?cos(2/C,  +  0,)] 

-- [f  i^2sCos(2A:i+2A:2s+<J>i)  -  5f  £fcos(2^i  +  0i)]  (3.23) 

gj  5SP 

—q^[EiE2E2SCOS(2Ki+K2+K2S+<Pi)-S^^EIEI^E2SCOS(2Ki-K2+K2S+<Px)] 

gj  cPS 

[^l£2£2SCOS(2£’i+^:2+^:2S+<i>l)-5f  £f£2^2sCOS(2^1+A:2-A'2S+<Pl)]  . 

To  form  the  Jacobian  matrix  for  this  pair  of  equations,  derivatives  of  FI  and  F2 
with  respect  to  k„  and  £n  are  needed.  These  derivatives  are  quite  complicated  due  to  the 
many  products  and  quotients  involving  the  it's  which  depend  on  both  kn  and  £n.  Also  the 
reflection  and  transmission  coefficients  depend  on  x's.  Since  the  inverse  of  the  Jacobian  is 
only  a  focusing  mechanism  to  the  next  guess,  the  derivatives  were  simplified  by  treating  the 
coefficients  as  constants.  The  penalty  for  this  is  a  probable  increased  number  of  iterations 
necessary  to  find  the  eigenvalues. 

Newton’s  method  works  well  if  the  slope  of  a  function  is  steq)  enough  near  a  root 
and  the  initial  guess  is  close  enough.  When  working  with  complicated  functions  such  as 
(3.22)  and  (3.23),  especially  when  it  is  not  easy  to  sketch  the  functions,  finding  the  roots 
becomes  a  difficult  task.  Initial  guesses  must  be  very  close  or  roots  will  be  missed. 

The  roots  found  by  the  factoring  technique  discussed  below  were  used  as  the  initial 
guesses.  Since  these  were  not  always  good  enough  to  find  all  the  exact  roots,  the  program 
also  used  guesses  in  between  the  factored  roots.  Due  to  slow  convergence,  round-off 
error,  and  insufficient  guesses,  some  roots  were  still  missed.  See  Appendix  for  program 
listing. 

To  determine  if  the  factored  eigenvalues  generated  by  the  program  were  reasonable, 
the  results  were  compared  to  those  produced  by  RAYMODE  using  the  same  input 
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parameters.  RAYMODE,  a  commly  used  propagation  loss  program,  was  altered  to  have 
the  same  reflecdon  coefficients  as  used  here  and  the  agreement  of  the  two  programs  was 
very  good. 

3.3.  Factoring  the  LEE  Characteristic  Equation 

The  factoring  technique  described  in  section  1.3  can  be  used  on  (3.21)  with  the 
introduction  of  some  product  errors.  Resulting  are  five  factors  to  be  called  the  water  and 
sediment  pp,  ss,  sp,  and  ps  factors. 

3.3.1.  The  LEE  Water  Equation 

From  (3.21),  the  water  equation  for  the  LEE  case  is 

.  (3.24) 

Taking  the  real  and  imaginary  parts  to  form  two  new  equations  leads  to  the  water  solutions. 


1 

,  , 

II 

(3.25) 

(3.26) 

The  T  -FF.  water  solutions  differ  firom  the  LLL  water  solutions  only  in  the  factor. 

3.3.2.  The  LEE  Sediment  Equations 

From  (3.21),  the  pp,  ss,  sp,  and  ps  factors  can  be  written.  The  method  for  solving 
each  equation  is  the  same.  Take  the  real  and  imaginary  parts  of  the  equation  to  form  two 


equations.  Solve  the  imaginaiy  equation  for  kn  using  Newton’s  method  for  one  variable 
and  substitute  the  result  into  the  real  equation  solved  for  £„.  The  solutions  to  the  sediment 
factored  equations  are 


-Inn  =0  ,  (3.27) 


1/2 


ln(5RiSPP) 


(3.28) 


F^Hkn)  =  lh{k\-kl)^^+  2r(it|s-it^)^/2+0j  _  2/i;r  =  0  ,  (3.29) 


ln(3liSSS) 


(3.30) 


F^^kn)  =  F^\kn)  *  2/i(/:f-/:2)‘/2+r(/:|-/t2)V2^.,(^2^^2)V2^.  0^_2nn=O,  (3.31) 


ln(9liSSP) 


J _ 


(3.32) 


— 


In  (91 15  PS) 


+/k^) 


(3.33) 


3.4.  Comparison  of  Factored  to  Exact  Solutions 

As  with  the  LLL  model,  extreme  parameter  sets  were  chosen  to  force  the  exact  solutions  to 
either  the  water  or  various  sediment  path  solutions.  For  the  force  to  water  test,  a  high 
density  and  sound  speed  are  used.  Figure  3.7  shows  the  agreement  between  the  water 
solutions  above  the  critical  angle  and  the  exact  solutions.  To  encourage  the  sedimrat 
paths,  similar  sound  speeds  for  water  and  sediment  are  entered  for  all  three  runs.  For  the 
pp  path,  a  small  sediment  attenuation  and  a  hi^  sediment  shear  attenuation  result  in  the 
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exact  solutions  being  equivalent  to  the  pp  solutions  as  shown  in  Figure  3.8.  Increasing  the 
pressure  wave  attenuation  for  the  sediment  while  decreasing  the  shear  wave  attenuation 
enhanced  the  shear  paths  as  shown  in  Figme  3.9. 


LEE  Forcing  Water  Solutions  to  Exact  Solutions 


■  water 
o  sediment  pp 

♦  sediment  ss 

•  sediment  ps/sp 

X  exact 


0.00  0.02  0.04  0.06  0.08  0.10  0.12  0.14 

Re(eigen  value) 


Hgure  3.7.  LEE  Force  of  Water  Solutions  to  Exact  Solutitxis 
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LEE  Forcing  Sediment  PP  Solutions  to  Exact  Solutions 


■  water 

o  sediment  PP 

■4-  sediment  ss 

•  sediment  ps/sp 

X  exact 


0.00  0.02  0.04  0.06  0.08  0.10  0.12  0.14 

Re(eigen  value) 


Figure  3.8.  LEE  Force  of  Sediment  PP  solutions  to  Exact  Solutions 


LEE  Forcing  Shear  Solutitms  to  Exact  Stdutions 


■  water 
o  sediment  PP 

*  sediment  ss 

•  sediment  ps/sp 

X  exact 


0.00  0,02  0.04  0.06  0.08  0.10  0.12  0.14 

Re(eigenvaiue) 


Hgure  3.9.  LEE  Force  of  Sedimoit  Shear  Solutions  to  Exact  Solutions 


Chapter  4 


RESULTS 


A  comparison  of  the  eigenvalues  generated  by  the  exact  solutions  and  the  factored 
soludons  in  both  the  LLL  and  T.FF  models  was  conducted  for  four  types  of  sediments:  a 
fine  fluid,  a  medium,  a  hard,  and  a  semiconsolidated  sediment  Some  typical  sediments 
that  fall  into  these  classes  are  silty  clay,  dll  and  sand,  coarse  sand  and  gravel,  and  chalk  and 
shale.  Table  4.1  lists  the  numbers  used  for  the  sound  speeds  and  attenuadon  coefficients 
for  the  water,  sediments  and  bedrock  (Hamilton  1980,  pp.  1326)  (Stoll  1986,  pp.  429) 
(Clay  and  Medwin  1977,  pp.  258)  (Beebe  1981,  pp.  78, 94, 1 15,  and  130). 

Two  criteria  are  used  to  evaluate  the  facuned  solution  eigenvalues:  placement  on  the 
eigenvalue  plot  of  -log(Iml£;,l)  vs.  kn ,  and  a  mode  summation  comparison  between 
factored  solutions  and  exact  solutions.  An  approximate  pressure  mode  summation  equation 
(4.1),  adapted  from  Clay  and  Medwin’s  (9.2.29),  is  an  incoherent  sum  for  a  source  and 
receiver  located  at  the  same  depth.  The  depth,  z,  is  100m  and  the  horizontal  range,  r,  is 
10km  fcH*  the  calculations.  A  log  ratio  of  the  water  (w)  plus  sediment  pp  (s)  solutions  sums 
divided  by  the  exact  (e)  solutions  multiplied  by  10  gives  the  nx>de  sum  comparison  in  dB 
as  in  (4.2).  The  closer  the  nxxle  sum  is  to  zero,  die  better  matched  the  factored  solutions 
are  to  the  exact  solutions. 

„  lcos(icS^^z)sitt(KS^^z)e-l^ 
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Table  4.1. 

Geoacoustical  Parameters  for  Water,  Sediments,  and  Bedrock 

Layer 

P 

c 

a 

c% 

as 

(g/cm3) 

(m/s) 

(dB/m  kHz) 

(m/s) 

(dB/m  kHz) 

Water 

1.0 

1500 

Fine  Fluid 

1.4 

1470 

0.01 

250 

15.0 

Medium 

1.9 

1600 

0.03 

450 

13.0 

Hard 

2.0 

1720 

0.003 

650 

10.0 

Semiconsolidated 

2.1 

2400 

0.05 

1000 

1.0 

Bedrock 

Q 

5500 

0.05 

2900 

(0.07) 

mode  sum  comparison  =  10  log  •  (4-2) 

Me 

4.1.  The  LLL  and  LEE  Reflection  and  Transmissions  Coefficients 

It  is  useful  to  examine  the  plots  of  die  reflection  and  transmission  coefficients  in 
Figures  4. 1-4.5,  made  for  a  medium  sediment  at  30  Hz.  Some  of  the  coefficients  affect  the 
eigenvalues  drastically  so  it  is  important  to  see  their  behavior. 

Figures  4. 1  and  4.2  compare  the  LLL  9ti2  and  9123  ^th  the  LEE  SRi2^^  and 
9t23^^-  There  is  litde  difference  between  the  two  Sin’s,  the  highest  difference  being  for 
the  lowest  modes,  or  highest  ’s.  SI23  is  equal  to  one  over  most  of  the  kn  range 
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R12 1  vs.  Re(eigenvalue) 


IR12I 

IR12ppl 


Figure  4. 1 .  Comparison  of  9li2  and  LEE  SR  12**^ 


I R23 1  vs.  Re(eigenvalue) 


IR23I 

IR23ppl 


Figure  4.2.  Comparison  of  LLL  9t23  and  LEE  9^23^ 


while  is  not  constant  or  even  smooth.  At  the  lowest  modes  the  falloff  of  9123*’*’  is 
most  important.  The  swings  in  5R23****  affect  the  shape  of  the  LEE  sediment  pp  solutions 
more  than  the  exact  solutions. 
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Figures  4.3-4.5  display  the  rest  of  the  LEE  coefficients.  Figure  4.3  shows  the 
extreme  difference  between  9i2i^*’  and  9l2i^^  for  the  LEE  case.  From  this  plot,  it  is  easy 
to  see  how  the  shear  waves  are  attenuated  quickly.  Figure  4.4  shows  the  LEE  9^23^^, 
9^23^*’  and  9^23*’^  coefficients.  They,  like  SR23*’*’,  are  not  smooth,  but  their  peaks  and  dips 
fall  in  the  region  of  low  kn .  Figure  4.5  displays  the  LEE  transmission  coefficients. 


Figure  4.3.  Comparison  of  LEE  9I21*’**  and  9l2i^^ 


Tiansmissions  Coefficients  vs.  Re(eigenvalue) 


IT12ppl 

IT12psl 

IT21ppl 

IT21spl 


Figure  4.5.  Comparison  of  LEE  Siz****.  3 12*^.  32i****  and  321*** 
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4.2.  LLL  Model  Results 

The  first  study  compares  not  only  the  differences  between  the  four  types  of 
sediments,  but  also  tests  variations  with  frequency  and  sediment  thickness.  Four 
combinations  of  thickness  and  frequency  are  examined:  low  thickness  (10m)  and  low 
frequency  (lOHz),  low  thickness  and  high  frequency  (50Hz),  high  thickness  (200m)  and 
low  frequency ,  high  thickness  and  high  frequency. 

4.2.1.  General  Trends 

From  (2.35),  the  expression  for  the  number  of  modes  in  a  duct,  it  is  evident  that 
with  increased  frequency  or  water  depth  there  will  be  an  increased  number  of  modes. 
Consequently  frequencies  above  50  Hz  were  not  studied  to  keep  the  number  of  modes 
manageable.  The  water  depth  for  this  study  was  arbitrarily  kept  constant  at  IQOtn. 
Comparing  Figures  4.6  and  4.7  with  4.8  and  4.9  shows  the  expected  increase  in  mode 
number  for  a  frequency  increase. 

For  the  larger  sediment  thickness ,  the  sediment  mode  attenuation  is  greater  than  for 
the  smaller  thickness  case.  The  wave  travels  further  in  the  sediment  and  is  therefore  more 
attenuated.  Also,  there  is  more  attenuation  in  the  higher  frequency  runs  because  the 
attenuation  factor  A  goes  up  with  frequency. 

4.2.2.  Comparison  of  Factored  Solutions  to  Exact  Solutions 

As  can  be  seen  in  Figures  4.6-4.9,  the  mode  attenuation,  or  1^1,  is  greater  for  both  factored 
solutions  than  for  the  exact  solutions.  One  of  the  factored  solutions  should  have  nearly  the 


Medium 
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Fine  Fluid  Medium 
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(|an|eAua3ia|  uii)3o|-  (lanieAuaSiaj  uii)3o|- 


(laniuAuaSpI  uii)8o|*  (|an|8Aua8{a|  uii)8o|> 


0.00  0.05  0.10  0.15  0.20  0.25  0.00  0.05  0.10  0.15  0.20  0.25 

Re(eigenvalue)  Re(eigenvalue) 

Figure  4.7.  LLL  Water,  Sediment,  New  Sediment,  and  Exact  Solutions  for  Four  Sediments:  t=10m,  f=50Hz 


Fine  Fluid  Medium 
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(laniUAuaSidl  uii)So|-  (|9n|eAU93i9|  uii)3o|> 


(laniBAuaSiaj  uii)3o|-  (laniBAuaSpI  aii)3oh 


Re(eigenvalue)  Re(eigenvaiue) 

Figure  4.8.  t  n  Water,  Sediment,  New  Sediment  and  Exact  Solutions  for  Four  Sediments:  t=200m,  f=I0Hz 


Fine  Fluid  Medium 
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(|9n|eAua3!9|  uii)3o|> 
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Re(eigenvalue)  Re(eigenvalue) 

Figure  4.9.  LLL  Water,  Sediment,  New  Sediment,  and  Exact  Solutions  for  Four  Sediments:  t=200m,  f=50Hz 
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equivalent  mode  attenuation  of  the  exact  eigenvalues  so  it  may  be  used  in  place  of  the  exact 
solutions.  From  Table  4.2,  the  mode  sum  comparisons  are  too  high,  as  high  as  30  dB  for 
the  semiconsolidated  lOHz,  10m  case. 

Since  the  water  factor  yields  solutions  that  are  straight  forward  to  calculate,  it  will 
not  be  changed.  An  adjusted  sediment  factor  is  needed. 

4.2.2. 1  Derivation  of  the  New  LLL  Sediment  Factor 

The  sources  of  error  in  the  sediment  factor  are  the  product  term  added  to  the  left 
side  of  the  characteristic  equation  as  described  in  section  1.4  and  setting  S2  equal  to  zero  in 
section  2.3.2.  Keeping  S2  in  the  equation  did  not  improve  the  results.  To  change  the 
sediment  factor,  keep  the  water  factor,  and  yet  remove  the  error,  a  new  sediment  factor  was 
developed  by  dividing  the  characteristic  equaticm  by  the  water  facton 

(1  - 

-  .  (4.3) 

(1 -9li9li2e-2^«) 

Writing  the  denominator  as  a  geometric  series,  expanding  that  series,  and  multiplying  gives 
(1  -  91i5li3e-2‘-^>  +  (4.4) 

After  substitution  for  9li3,  the  positive  terms  in  9li2  drop  out  leaving 

1 - [1  +  +  . .  .]^-2i(iir,  (4.5) 

1  —  9l2l9l23A^e"2iifi 

The  new  sediment  factor  is  to  be  of  the  form 

1-C«-2^^>  +  *^2),  (4.6) 

where  C,  with  the  denominator  expanded,  is 
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C  =  ^iSn^23^2iA^  [l  +  SRi9li2e-2**^i  +  +•••]>< 

[l  —  +  SR  1 29^23^ - ]  • 

In  order  to  obtain  a  smooth  mode  attenuation,  2Ari  is  approximated  as  much  larger  than 
2K2,  which  amounts  to  h  being  much  larger  than  t.  The  eigenvalues  of  this  factor  with 
combined  phase  ^ri+^r2  will  be  very  close  to  those  of  the  water  factor,  since  the  sediment 
phase  is  very  small.  Both  exponentials  and  e-2‘^2  are  approximated  as  unity  in  the 
magnitude  of  C. 

|Cl  =  SRi3i29l2332i>i2[l  +SRiSRi2  +  Sk5sr52][i  -9^129^23/^2 +  9t529?23A^]  .  (4.8) 

Both  series  work  best  truncated  to  quadratic  form  as  in  (4.  8).  Higher  orders  gave  mode 
attenuations  of  the  wrong  sign,  indicating  the  sums  were  diverging.  The  original  sediment 
factor  has  the  two  terms  in  brackets  equal  to  unity. 

4. 2. 2. 2  Comparison  of  New  LLL  Sediment  Factor  with  Exact  Solutions 

From  Figures  4.6-4.9  the  new  sediment  factor  ^’s  arc  much  closer  to  the  exact 
£V,’s.  Table  4.2  confirms  the  improvement  with  much  smaller  mode  sum  comparisons. 

The  new  sediment  factor  works  very  well  for  the  fine  fluid,  medium,  and  hard  sediments, 
but  is  not  improved  enough  for  the  semiconsolidated  sediment 

4.3.  LEE  Model  Results 

The  tests  cm  the  l^FE  model  were  done  for  a  constant  low  frequency  (30Hz)  since 
the  affects  of  frequency  have  already  been  shown  in  the  LLL  study.  Two  sediment  depths, 
lOm  and  SOOm,  were  examined  for  a  water  depth  9Q0m.  The  results  are  shown  in  Figures 
4.10-4.11  and  Table  4.3. 
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Table  4.2. 


LLL  Mode  Sum  Comparisons  in  dB  with  Sediment  and  New  Sediment  Factors 


Fine  Fluid 


Medium 


Hard 


Semiconsolidated 


Fine  Fluid 


Medium 


Hard 


Semiconsolidated 


10m,  lOHz 

10m,  50Hz 

Sediment  New  Sediment 

Sediment 

New  Sediment 

-2  0 

-2 

0 

-13  -1 

-11 

0 

-21  -2 

-19 

-2 

-30  -9 

-14 

-9 

200m.  lOHz 

200m,  50Hz 

Sediment  New  Sediment 

Sediment 

New  Sediment 

-2  0 

-1 

1 

-9  -1 

-7 

1 

-15  -1 

-14 

0 

-27  -9 

-A 

-4 
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There  are  four  sediment  factors,  pp,  ss,  sp,  and  ps.  The  pp  mode  attenuation  has  a 
dip  and  peak  corresponding  to  the  dip  and  peak  in  the  9^23^^  plot  in  Figure  4.2.  The  ss, 
ps,  and  sp  mode  attenuations  are  much  lower  than  the  pp.  The  ps  and  sp  solutions  are 
actually  equivalent  since  their  coefficient  products,  3i2**®9l23***32i  and  3i29l23*’^32i^’’ 
are  equivalent. 

The  water  and  sediment  pp  solutions  are  the  best  approximation  to  the  exact 
solutions.  The  ss  and  ps/sp  solutions  are  usually  highly  attenuated  since  3 12**^  and  321^** 
are  small,  and  can  be  left  out  of  the  LFF  model.  The  dips  and  peaks  in  the  sediment  pp 
solutions  are  cause  for  concern  when  the  using  the  factors  in  place  of  the  exact. 

The  program  missed  some  of  the  modes  for  the  LEE  exact  case.  The  accuracy  of 
mode  sum  con^arisons,  which  include  only  the  water,  sediment  pp,  and  exact  solution 
sums  as  in  the  LLL  case,  is  thus  somewhat  limited. 

A  new  LEE  sediment  pp  was  developed  analogous  to  the  one  for  the  LLL  new 
sediment  factor.  With  both  series  linear,  quadratic,  cubic,  fourth  order  and  various 
combinations  of  those,  £^’s  of  the  wrong  sign  wa?e  generated.  Also  the  mode  attenuation 
plot  was  not  smooth,  so  the  method  of  coiTection  was  abandoned. 

4.4.  Comparison  of  LEE  Model  to  LLL  Model 

Figures  4.12-4.13  compare  the  i  -U- and  LFR  water  and  sediment  pp  solutions. 
The  water  factors  for  the  LEE  and  LU.  models  are  nearly  identical  for  identical  input 
parameters.  This  is  expected  since  is  the  same  fn*  both  and  for  the  liquid/elastic 
boundary  is  nearly  equal  to  91 12  for  the  liquid/liquid  boundary  as  shown  in  Figure  4. 1 . 
Consequently,  the  LLL  water  factor  can  be  in  a  LEE  model. 


Fine  Fluid  Medium 
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(laniUAuaSial  uii)8o|-  (|an|«Aua8|a|  uii)8o|- 


Re(eigenvalue)  Re(eigenvalue) 

Figure  4. 10.  LEE  Water,  Sediment  PP,  SS,  and  PS/SP,  and  Exact  Solutions  for  Four  Sediments:  t=10m,  f=30Hz 


Fine  Fluid  Medium 
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(IsnieAUsSial  iui)So|*  (lanieAUdSfal  uii)3o|* 


(|an|BAua8ia|  uii)8o|*  (|an|BAua8|a|  uii)8o|- 


Re(eigenvalue)  Re(eigenvalue) 

Figure  4.1 1.  LEE  Water,  Sediment  PP,  SS,  and  PS/SP,  and  Exact  Solutions  for  Four  Sediments:  t=300m,  f=30Hz 
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Table  4.3. 

Approximate  T  .F.F.  Mode  Sum  Comparisons  in  dB 

10m,  30Hz 

300m,  30Hz 

Fine  Fluid 

6 

2 

Medium 

8 

10 

Hard 

8 

3 

Semiconsolidated 

16 

7 

The  sediment  solutions  do  not  propagate  as  well  in  the  LEE  environment  as  they  did 
in  the  LLL  environment  due  to  the  added  shear  loss.  9^23^  not  equal  to  one  for  the  LEE 
case  as  it  was  for  the  LLL  case  as  shown  in  Figure  4.2,  but  instead  has  dips  and  peaks, 
thus  changing  the  shape  of  the  sediment  pp  solutions. 


Medium 
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Chapter  5 

CONCLUSIONS 

Approximate  solutions  for  the  three  layer  duct  characteristic  equation  were 
developed  for  quickly  obtaining  the  eigenvalues.  The  LLL  model  factors  are  both  quick 
and  a  good  approximation  of  the  T.T.T.  exact  solutions.  The  T.FF.  factored  solutions  are  also 
easy  to  calculate,  but  the  loss  of  accuracy  may  outweigh  the  benefits  of  speed. 

The  method  of  factoring  described  in  section  1.4  was  applied  to  both  the  LLL  and 
LEE  characteristic  equations.  The  factored  solutions  were  then  compared  to  the  exact 
solutions  for  each  model.  Also  the  l.i.L  and  LEF.  factored  solutions  were  compared  to  each 
other.  Variations  in  sediment,  sediment  thickness,  and  frequency  were  examined. 

The  LLL  faaored  solutions  predicted  losses  much  higher  than  the  LLL  exacts.  A 
new  sediment  factor  was  then  develqjed,  based  cm  the  characteristic  equation  and  the  water 
factor.  The  losses  for  the  new  factor  are  much  closer  to  the  losses  of  the  exact  solutions 
with  the  best  results  for  the  fine  fluid,  medium,  and  hard  sediments.  The  losses  of  the  new 
factcH*  for  the  semiconsolidated  sediment,  which  favors  the  water  path,  were  improved,  but 
not  as  much  as  the  others. 

The  factored  solutions  for  the  LEE  nnodel  also  predicted  mcne  losses  than  the  LEE 
exact  solutions.  Attempts  to  correct  the  sediment  pp  factor  proved  inadequate.  Even 
though  this  I  FF  factoring  scheme  predicts  too  much  loss,  using  it  may  still  be  better  than 
using  a  T.T.I.  model  which  would  predict  far  too  little  loss  for  an  elastic  sediment 
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APPENDIX 

COMPUTER  PROGRAM  LISTING 


The  following  computer  listing  is  a  program  that  finds  the  eigenvalues  for  the  LLL 
water,  sediment,  new  sediment,  and  exact  solutions.  Also  it  finds  the  LEE  water,  sediment 
pp,  ss,  sp,  and  ps,  and  exact  solutions.  The  program  is  written  in  HP  BASIC.  The 
symbols  and  variable  names  used  in  the  program  do  not  necessarily  match  the  symbols 
used  in  the  text. 


RAD 

OPTION  BASE  1 
GINIT 
GCLEAR 
GRAPHICS  ON 

DIM  Kw(90),  Ew(90),  Ks(90),  Es(90),  Ksm(90),  Esm(90),  Kn(90),  En(90) 
DIM  Kw  1  (90),  Ew  1  (90),  Kw2(90),  Ew2(90) 

DIM  Kw3(90).  Ew3(90),  Kw4(90),  Ew4(90),  Knl(90),  Enl(90) 

COM  /Parameterl/  Cl,  C2,  Rhol,  Rho2,  H,  T,  Freq,  Alfa2,  Rl,  Rlphase 

COM  /Parameter2/  C2s,  C3,  C3s,  Rho3,  Alfa2s 
COM  /Wavenumber/  W,  Kl,  K2,  K2s,  K3s 

Cl=1500 

C3=2900 

Rhol=1.0 

Rho3=2.7 

H=900 

T=10 

Freq=30 

Rl=0.99 

Rlphase=PI 

CALL  Label4axes(l) 

F0RDiffcase=lT04 
SELECT  Diffcase 
CASEl 

C2=1470 

C2s=250 

Rho2=1.4 
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Alfa2=0.01 
Alfa2s=15.0 
CASE  2 

C2=1600 
C2s=450 
Rho2=1.9 
Alfa2=0.03 
Alfa2s=13.0 
CASE  3 

C2=1720 
C2s=650 
Rho2=2.0 
Alfa2=0.03 
Alfa2s=10.0 
CASE  4 

C2=2400 
C2s=1000 
Rho2=2.1 
Alfa2=0.05 
Alfa2s=1.0 
END  SELECT 


W=2*PI*Freq 

K1=W/C1 

K2=W/C2 

K3sW/C3 

K2s=W/C2s 

K3s=W/C3s 

ALPHA  ON 
OUTPUT  KBD;“  K”; 


CALL  Water(Kw(*),Ew(*),Nmodew“IirV‘o”,Nofextras) 

CALL  Watcr(Kw  1  (*),Ew  1  (*)  J4mode  1  w,“lee’V‘o”,Nofextras) 

CALL  Scd(Ks(*),Es(*),Kw(*),Nmodes,Nmodcw,“in”,“  pp”  “+”,Nofcxtras,“s  ”) 
CALLSed(Ksm(*),Esm(*),Kw(*),Nmodcms,Nmodew,“llm”.“pp”,“u”,Nofcxtras,“sm”) 
CALL  Scd(Ks l(*),Es l(*),Kw l(*),Nniodcls,Nmodel w,“lee”,“  pp”,“+”,Nofcxtras,“sl ”) 
CALLSed(Ks2(*),Es2(*),Kwl(*),Nmodc2s,Nmodclw,“lcc”.“  ss”,“*”,Nofextras,“s2”) 
CALLSed(Ks3(*),Es3(*),Kwl(*),Nmodc3s,Nmodclw,“lec”,“sp”"x”,Nofextras“s3”) 
CALL  Sed(Ks4(*),Es4(*),Kw  1  (*),Nmodc4s,Ninodc  1  w,“lec”,‘‘  ps”,“-”,Nofextras,“s4”) 

CALL  Exaci(Kn(*)£n(*),Kw(*)^w(*),Ksin(*),Esin(*),Ninoden,Nmodcw,Nmodcms, 
»lll..«A»NofextrasJ5ifFcase) 

CALL  Exact(Kn  1  (*),En  1  (*),Kw  1  (*),Ew  1  (*),Ksl  (*)^s  1  (*).Nmodel  n.Nmodcl  w, 
Nmcxlc  1  s,“lce”,“n”,Nofcxtras4>iffcasc) 

CALL  Sort(Kn('*')£n(*),Naiodcn) 

CALL  Sort(Knl(*)^nl(*),Ninodeln) 

CALL  Nat_log_c(Ew(*)^n_cw(*)^tnodcw) 

CALLNat_log_e(Es(*)»Ln_cs(*),Nmcxlcs) 

CALLNatJog_c(Esm(*)^n_csin(*),Ninodcnis) 
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CALL  Nat_log_e(En(*),Ln_en(*),Nniodcn) 

CALL  Nat_log_e(Ew  1  (*),Ln_ew  1  (*),Nmocic  1  w) 

CALL  Nat_log_e(Es  1  (*)  Xn_es  1  (*  ),Nmode  1  s) 

CALL  NatJog_e(Es2(*)4-n_es2(*),Nmode2s) 

CALL  Nat_log_e(Es3(*)J-n_es3(*),Nmode3s) 

CALL  Nat_log_e(Es4(*)4-n_es4(*),Nmodc4s) 

CALL  Nai_log_e(En  1  (*),Ln_en  1  (*),Nmodel  n) 

CALL  Fouraxes(5  JJiffcase) 

CALL  Ln_e_vs_k_plot(Kw(*),Ew(*),Ln_ew(*),Nmodew,5,”o”) 

CALL  Ln_e_vs_k_plot(Ks(*)£s(*)4-n_es(*),Nmodes,l 

CALL  Ln_e_vs_k_plot(Ksni(*),Esm(*),Ln_esm(*),Nmodems,4,”u”) 

CALL  Ln_e_vs_k_plot(Kn(*)^n(*),Ln_en(*),Nmoden,7,”'^”) 

CALL  Ln_c_vs_k_plot(Kw  1  (*),Ew  1  (*),Ln_cw  1  (♦),Nmo(le  1  w,6,”o”) 

CALL  Ln_e_vs_k_plot(Ks  1  (*)  JEs  1  (*)a.n_esl(*),Nmode  1  s,3,”+”) 
CALLLn_e_vs_k_plot(Ks2(*),Es2(*),Ln_es2(*).Nmode2s,4”*”) 
CALLLn_e_vs_k_plot(Ks3(*).Es3(*)Xn_es3(*).Nm(xie3s,5,”x”) 
CALLLn_e_vs_k_plot(Ks4(*),Es4(*).Ln_es4(*),Nmode4s,7”-”) 

CALL  Ln_e_vs_k_plot(Kn  1  (*),En  1 (*),Nmode  1  n,2,”n”) 

PRINT  “LLL  Mode  Sum” 

CALL  Incohmt(Kw(*),Ew(*),Ks(*),Es(*),Kn(*),En(*),Nmodew,Nmodes,Nmoden, 
Indbex.Nofextras) 

PRINT  “LLL  Mode  Sum  with  New  Sediment  Factor” 

CALL  Incohmt(Kw(*),Ew(*),Ksm(*),Esm(*),Kn(*),En(*),Nmodew,Nmodems, 
Nmoden.Indbex.Nofextras) 

PRINT  “LEE  Mode  Sum” 

CALL  Incohmt(Kw  1  (*),E w  1  (*),Ks  1  (♦),Es  1  (*),Kn  1  (*),En  1  (*),Nmodel  w.Nmode  1  s, 
Nmode  1  n, Indbex.Nofextras) 

NEXTDiffcase 

END 


SUB  Water(Kw(*),Ew(*),Nmodew,Type$,Sym$,Nofextras) 

CX»4  /Parameterl/  Cl,  C2,  Rhol,  Rho2,  H,  T,  Freq,  Alfa2,  Rl,  Rlphasc 
COM  /Wavenumber/  W,  Kl,  K2,  K2s,  K3s 

Frmtl:IMAGE4X,4A,2D,4A,SD.4DE,5X,4A,2D,4A,SD.4DE,3X.D.2D,5X,2D.D 
INTEGER  N,M,I 

FOR  1=1  TO  90 
Kw(I)=0.0 
NEXT  I 
PRINT 

PRINT  TAB(  19);“WATER  ”;Type$;“  SOLNS”;Sym$;“  R  Phil ” 

PRINT 

Nofextras=0 


63 

Nmodew=0 

N=1 

WHILE  ((2*N*PI-Rlphase)/(2*H))<=K1 

Kw(N)=SQR(Kl''2-((((2*N*PI-RIphase)/(2*H))''2) 

Kappa  1  =SQR(K  1  '^2-Kw(N)^2) 

IF  (K2^2-Kw(N)'^2)<0  THEN 
Kappa2=0 

Noextras=Noextras+ 1 
ELSE 

Kappa2=SQR(K2'^2-Kw(N)'^2) 

END  IF 

Phil=ACS(Kw(N)/Kl)*  180/PI 
Z12=(Rhol  *Kappa2)/(Rho2*Kappal) 

R1211=(l-Z12)/(1+Z12) 


Nmodew=Nmodew+ 1 

M=Nmodew 

Kw(M)=Kw(N) 

SELECT  Types 
CASE  **111” 

Ew(M)=LCX3(ABS(R  1  *R  1 211))*Kappa  l/(Kw(M)*2*H) 

PRINT  USING  Fraitir  Kw(”,M,“)  =  ",Kw(M),“  Ew(”,M,“)  «  ”^w(M). 
RI211.Phil 

CASE  “lee” 

IF  (K2s''2-Kw(N)'^2)<=0  THEN 
R121e=l 
ELSE 

Kappa2s=SQR(K2s''2-Kw(N)''2) 

Phi2s=ATN(Kappa2s/Kw(]^)*  1 80/PI 

R12num=(-Z12+l-(SIN(2*Phi2s))''2*(l-Kappa2/Kappa2s)) 

R12den=(Z12+l-(SIN(2*Phi2s)y'2*(l-Kappa2/Kappa2s)) 

R 1 21esR  1 2nuin/R  1 2den 
END  IF 

Ew(M)=LOG(ABS(Rl*R121c))*Kappal/(Kw(M)*2*H) 

PRINT  USING  Frmtl;“Kwl(”,M,“)  =  ”JCw(M),“EwlC’,M ,“)  =  ”£w(M), 
R121e,Phil 

END  SELECT 
N=N+1 
END  WHILE 
PRINT 
SUBEND 


SUB  Scd(Ks(*)3s(*),Kw(*),Nmodes,Nmodew,Type$,Kind$,Sym$,Extras,S$) 
Frmtl:  IMAGE  4X,SD.3DE,6X.2D,SD.3DE,5X,D.3D,3X,D.2D,3XX>.2D.4X.2D.D 
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FnTit2:  IMAGE  4X,SD.3DE,6X,2D,SD.3DE,5X,D.3D,2X,2D.2D,3X,2D.2D,3X, 
2D.2D.3X.2D.D 

FrmtS:  IMAGE  8X,A,2A.10X,A,10X,A.2A,9X,5A,3X,3A.4X.3A,4X,4A 
Frmt6:  IMAGE  8X,A,2A.10X,A,10X,A,2A.9X,5A,3X.3A,5X,3A,5X,3A,3X,4A 
Frmt?:  IMAGE  8X.A,2A,10X,A,10X,A,2A,9X,5A,2X.5A,2X,5A.2X,5A,3X.4A 

CX3M  /Parameterl/  C 1 .  C2.  Rho  1 .  Rho2.  H,  T.  Freq,  Alfa2.  R 1 ,  R Iphase 

COM  /Parameter2/  C2s,  C3,  C3s,  Rho3,  Alfa2s 

COM  /Wavenumber/  W,  Kl,  K2,  K2s,  K3s 

INTEGER  Mm,  Num_divs,  N,  M,  Stepslow,  Nn 

PRINT 

PRINT  TAB(16);“SEDIMENT  ”;Type$;“  SOLNS”;Sym$ 

PRINT 

SELECT  TYPES 
CASE“m” 

PRINT  USING  Frmt5:“K”,S$,“N”,“E”,S$,“Atten’V‘R12”,“R23”, “Phil” 

CASE  “Urn” 

PRINT  USING  Frmt5;“K”,S$,“N”,“E”,S$.“Atten”,“R12”,“R23”,“Phil” 

CASE  “lee” 

SELECT  Kinds 
CASE“  pp” 

PRINT  USING  Frmt6;“K”,SS,“N”,“E”,SS,“Atten”,“R23”,“T12pp”,“T21pp”,“Phi  1” 
CASE  “  ss” 

PRINT  USING  Frmt7;“K”,SS.“N”,“E”,SS.“Atten”,“R23ss”.“T12ps”,“T21sp”,“Phil” 
CASE“sp” 

PRINT  USING  Frmt7;“K”,SS,“N”,“E”,SS,“Atten”,“R23sp”,“T12ps”,“T21pp”,“Phil” 
CASE  “  ps” 

PRINT  USING  Frmt7;“K”,SS,“N”.“E”.SS,“Atten”,“R23ps”,‘T12pp”,“T21sp”, “Phil” 
END  SELECT 
END  SELECT 


Onum_divs=Nmodew-Extras 

Max_iters=40 

Nmodes=0 

FOR  1=1  TO  90 
Ks(I)=0.0 
NEXT  I 


Nn=l 

REPEAT 

X=MIN(0.9999*K  1 ,0.99999*K2) 

Kappa  1  =SQR(K  1  '^2-X'^2) 
Kappa2=SQR(K2''2-X''2) 
Kappa2s=SQR(K2s''2-X''2) 
Phase=2*Nn*PI-R  Iphase 

IFKindS=“pp”THEN 

Fs=2*H*Kappa  1  +2*T*Kappa2-Phase 


ELSE 


IF  Kind$=“  ss”  THEN 

Fs=2*H*Kappal+2*T*Kappa2s-Phase 

ELSE 

Fs=2*H*Kappal+T*(Kappa2+Kappa2s)-Phase 
END  IF 
END  IF 

Nn=Nn+l 

UNTIL  Fs<0  OR  Nn>=90 

Nstart=Nn-l 

N=Nstart 

Stepslow=I 

Rootsave=10000 

IF  Kind$=“  pp’THEN 

Tst=2*H*Kl+2*T*K2 

ELSE 

IF  Kind$=“  ss”  THEN 

Tst=2*H*Kl+2*T*K2s 

ELSE 

Tst=2*H*K  1 +T*  (K2+K2S) 

END  IF 
END  IF 

FOR  Stepslow=Stepslow  TO  Nmodew+ 1-Extras 
ff  Stepslow=l  THEN 

Kstan=Kw(StepsIow+Extras) 

Ks  ipl=MIN(0.9999*Kl,0.99999*K2) 

ELSE 

IF  Stepslow<Nmodew+ 1-Extras  THEN 
Kstan=Kw(Stepslow+Extras) 

Ks_ip  1  =Kstart=K  w(Stepslow- 1  +Extras) 

ELSE 

Kstart=l.E-4 

Ks_ip  1  =Kw(Stepslow- 1  +Extras) 

END  IF 
END  IF 

Num_divs=INT(Onum_divs*(  1  +Stcpslow/(Nmodew+l  ))/Freq*0. 1 )) 
Kstcp=(Ks_ip  1-Kstart)/Num_divs 

FOR  L=1  TO  Nutn_divs 
Num_iters=0 
Phase=2*  N*  PI-R 1  phase 
If  Tst<Phase  THEN  Sub_sed_end 
Ks_ip  1  =Kstan+(L- 1  )*Kstcp 

REPEAT 

Ks  i=Ks_ipl 

Kappal=SQR(Kl''2-Ks_i''2) 

Kappa2=SQR(K2'^2-KsJ''2) 

Kappa2s=SQR(K2s''2-KsJ''2) 


Signl=SGN(Kl''2-Ks_i'^2) 

Sign2=SGN(K2^2-Ks_i''2) 

Sign2s=SGN(K2s'^2-Ks_i'^2) 

IFKind$=“pp”THEN 

Fs=2*H*Kappa  1  *Sign  1  +2*T*Kappa2*Sign2-Phase 
Fsprime=-2*Ks_i*(H/(Kappal*Signl)+T/(Kappa2*Sign2)) 

ELSE 

IF  Kind$=“  ss”  THEN 

Fs=2*H*Kappal*Signl+2*T*Kappa2s*Sign2s-Phase 
Fspriine=-2*  Ks_i*  (H/(Kappa  1  *  S  ign  1  )+T/(Kappa2s*  Sign2s)) 

ELSE 

Fs=2*H*Kappal+T*(Kappa2’'‘Sign2+Kappa2s*Sign2s)-Phase 

Fsprime=-Ks_i*(2*H/(Kappal*Signl)+T/(Kappa2*Sign2) 

+T/(Kappa2s*Sign2s)) 

END  IF 
END  IF 

IF  ABS(Fs)>10000  THEN  Nxtl 

Ks_ip  1  =Ks_i-Fs/Fsprime 
IF  Ks_ipl>Rootsave  THEN  Ks_ipl=Rootsave 
N  utn_iters=N  um_iters+ 1 

UNTIL  ABS(Fs)<=LE-4  OR  Num_iiers>=Max_iters 

IF  Num_iters>=Max_iters  THEN 
PRINT  “Max  iterations” 

ELSE 

IF  KsJpl>MIN(0.9999*Kl,0.99999*K2)  THEN  Nxtl 

M=N-Nstart+1 

Ks(M)=KsJpl 

Rootsave=Ks_ipl 

FOR  Mm=l  TO  Nmodes 

Testy=ABS(Ks(M)-Ks(Mm)) 

IF  Testy<=0.00001  THEN  Nxtl 
NEXT  Mm 

Nmodes=Nmodcs+ 1 
IF  Nmodes=l  THEN  PRINT 

Phi  1  =ATN(Kappal/Ks(M)) 

Phi2=ATN(Kappa2/Ks(M)) 

Rng=T/SIN(Phi2) 

D=Rhol/Rho2 

PsKappa2/Kappal 

Z12=D*P 

Z23=(Rho2/Rho3)*(Kappa3/Kappa2) 

R12U*ABS((1-Z12)/(1+Z12)) 

R23U=ABS((1  -Z23)/(  1+Z23)) 


SELECT  Types 
CASE“U1” 


Atten=10^(-2*Rng*Freq*Alfa2/10000) 

S=Atten*(  1-R1211''2)*R2311 

PRINT  USING  Frmtl;Ks(M),M3s(M),Atten,R1211,R2311,Phil 
CASE  “11m” 

Anen=  10^(-2*Rng*Freq*  Alfa2/100(X)) 

S=Atten*(l-R1211^2)*R2311 

MagC=Rl*S*(l+Rl*R1211+RlA2*R1211'^2)*(l-R12U*R2311*Atten+ 

R1211''2*R2311''2*Atten''2) 

PRINT  USING  Frmtl;Ks(M),M^s(M),Atten.R1211,R2311,Phil 
CASE  “lee” 

SELECT  Kinds 
Phi3=ATN(Kappa3/Ks(M)) 

Phi2s=ATN(Kappa2s/Ks(M)) 

Phi3s=ATN(Kappa3s^s(M)) 

Rngs=T/SIN(Phi2s) 

B 1 2=  1  -(SIN(2*Phi2s))''2*(  1  -Kappa2/Kappa2s) 
R121e=(-Z12+B12)/(Z12+B12) 

T12=2*D*COS(2*Phi2s)/(Z12+B  12) 

T21=2*P*COS(2*Phi2s)/(Z12+B  12) 

T12ps=-^*D*(COS(Phi2s)'^2)*TAN(Phi2)/(Z12+B12) 

T21sp^*P*(COS(Phi2s)'^2)*TAN(Phi2s)/(Z12+B12) 

R21den=SIN(Phil)*((C2/C2s)*COS(2*Phi2s)'^2+SIN(2*Phi2)* 

SIN(2*Phi2s)*(C2s/C2))+(C2*Cl*Rhol*SIN(Phi2))/(C2s*C2*Rho2) 
R2 1  =((SIN(Phi  1  )*((C2/C2s)*COS(2*Phi2s)'^2-SIN(2*Phi2)* 
SIN(2*Phi2s)*(C2s/C2))-(C2*Cl  *Rhol  *SlN(Phi2))/ 
(C2s*C2*Rho2)))/R2Iden 

R2 1  ss=((SIN(Phi  1  )*((C2/C2s)*COS(2*Phi2s)'^2-SIN(2*Phi2)* 
SIN(2*Phi2s)*(C2s/C2))+(C2*Cl*Rhol*SIN(Phi2))/ 
(C2s*C2*Rho2)))/R21den 

R21ps=^2*SQR(SIN(2*Phi2)*SIN(2*Phi2s))*COS(2*Phi2s)*SIN(Phil)* 

SIN(2*Phi2)/SIN(2*Phi2s))*(C2s/C2)/R21den 

R21sp=^2*SQR(SIN(2*Phi2)*SIN(2*Phi2s))*COS(2*Phi2s)*SIN(Phil)* 

SIN(2*Phi2s)/SIN(2*Phi2))*(C2/C2s)/R21den 

Tan2=Kappa2/Ks(M) 

Tan3=Kappa3/Ks(M) 

Tan2s=Kappa2s/Ks(M) 

Tan3s=Kappa3s^s(M) 

Mu2s=Rho2*C2s''2 

Mu3s=Rho3*C3s''2 

Pl=(2*Mu2s+Mu3s*(Tan3s^2-l))*Tan2*Tan2s 

P2=(Mu2s*(Tan2s''2-l)-Mu3s*(Tan3s''2-l))*Tan2 

P3a:(2*Mu3s+Mu2s*(Tan2s''2-l))*Tan3*Tan2s 

P4=2*(Mu3&-Mu2s)*Tan3*Tan2*Tan2s 

Ql=2*(Mu2s-Mu3s)*Tan3*Tan2*Tan2s 

Q2=(2*Mu3s+Mu2s*(Tan2s''2-l))*Tan3s*Tan2 

Q3=((-Mu2s*(Tan2s''2-l))+(Mu3s*(Tan3s''2-l)))*Tan2s 

Q4«P1 
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Denomin=((P  1  +P3)*(Q2+Q4)-(P2+P4)*(Q  1 +Q3)) 
R23pp=(((Pl-P3)*(Q2+Q4)-(P2+P4)*(Ql-Q3))/Denomin) 
R23ss=(((P2-P4)*(Ql+Q3)-<Pl+P3)*(Q2-Q4))/Denomin) 
R23sp=(2*(P3*Ql-Pl  *Q3)/Denomin) 

R23ps=(2*(P4*Q2-P2*(^)/Denomin) 

CASE“pp” 

Coefs=Rl*R23*T12*T21 

Atten=  10^(-2*  Alfa2*Freq*Rng/l{X)00) 

Es(M)=LOG(ABS(Coefs*Atten))/(2*Ks(M)*(H/Kappal+T/Kappa2)) 
PRINT  USING  Fnnt2;Ks(M),M,Es(M),Atten,R23,T12,T21,Phil 
CASE  “  ss” 

Coefs=R  1  *R23ss*T  1 2ps*T2 1  sp 
Atten=I0^(-2*Alfa2s*Freq*Rngs/I00(X)) 

Es(M)=LOG(ABS(Coefs*Atten))/(2*Ks(M)*(H/KappaI+T/Kappa2s)) 
PRINT  USING  Fnnt2;Ks(M),M,Es(M),Atten,R23ss,TI2ps,T21sp,Phil 
CASE“sp” 

Coefs=R  1  *R23sp*T  1 2ps*T2 1 
Atten=10^(-(Alfa2*Rng+Alfa2s*Rngs)*Freq/100(X)) 
Es(M)=LC)G(ABS(Coefs*Atten))/(Ks(M)*(H/Kappal+T/Kappa2+ 
T/Kappa2s)) 

PRINT  USING  Fnnt2;Ks(M),M,Es(M),Atten,R23sp.T12ps,T21,Phil 
CASE“  ps” 

Coefs=Rl*R23ps*T12*T21sp 

Atten=10^(-(Alfa2*Rng+Alfa2s*Rngs)*Freq/10000) 

Es(M)=LOG(ABS(Coefs*Atten))/(Ks(M)*(H/Kappal+T/Kappa2+ 

T/Kappa2s)) 

PRINT  USING  Frmt2;Ks(M).M,Es(M),Atten,R23ps,T12,T21sp,Pbil 
END  SELECT 
END  SELECT 


N=N+1 
END  IF 

NxU;! 

NEXTL 

NEXT  Stepslow 
PRINT 
Sub_sed_end;! 
SUBEND 


SUB  Exact(Kn(*),En(*),Kw(*),Ew(*),Nino(len,Nmodew,Nmodes,Type$,Syni$,Extras, 

SS.OPTIONAL  INTEGER  Diffcase) 

Frmtl  :IMAGE  3X,4A,2D,4A,SD.4DE,4X,4A,2D,4A,SD.4DE,2X,2D.D 
Fnni2:IMAGE  3X,4A^D,4A,SD.4DE,4X,4A,2D,4A,SD.4DE,3X,2D.2D,2X,2D.2D, 
3X  3D  D 

Fnnt3:  IMAGE  2X,3D.2D.3X.2D.2D.3X,2D.2D,3X.3D.2D,2X.2D.2D,3X,2D.2D,3X, 
2D.2D,4X,2D.2D,2X,3D.2D 
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COM  /Parameierl/  Cl,  C2,  Rhol,  Rho2,  H,  T,  Freq,  Alfa2,  Rl,  Rlphase 

COM  /Parameter2/  C2s,  C3,  C3s,  Rho3,  Alfa2s 

COM  /Wavenumber/  W,  Kl,  K2,  K2s,  K3s 

INTEGER  I,  L,  M,  N,  Num_iters,  Max_iters,  Num_divs,  Stepslow 

DIM  Wi(2.1),  Wipl(2,l),  F(2,l),  1(2,2),  Jinv(2,2) 

PRINT 

PRINT  TAB(19);“EXACT  ”;Type$;“SOLNS”;Sym$ 

PRINT 

Max_uers=100 

Wipl(l,l)=0 

Wipl(2,l)=0 

Nmoden=0 

Stepslow=l 

Onum_divs=Nmodes/ 1 .2 

FOR  1=1  TO  90 
Kn(I)=0 
En(I)=0 
NEXT  I 

FOR  Stepslow=Stepslow  TO  Nmodes+1 
SELECT  Type$ 

CASE“m” 

IF  Stepslow=l  THEN 

Kstan=(Kw(Stepslow+Extras)+Ks(Stepslow+Extras))/2 

Estait=(Ew(Siepslow+Extras)+Es(Stepslow+Extras))/2 

Knext=MIN(0.9999*Kl,0.99999*K2) 

ELSE 

IF  Stepslow<Nmodes+l  THEN 

Kstart=(Kw(Stepslow+Extras)+Ks(Stepslow+Extras))/2 
Estan=(Ew(Stepslow+Exiras)+Es(Stepslow+Extras))/200 
Kncxt=Ks(Stepslow-l ) 

ELSE 

Kstan=l.E-4 

Estan=(Ew(Stepslow+Extras)+Es(Stepslow+Extras))/200 
Knext=Ks(Stepslow-l ) 

END  IF 
END  IF 
CASE  “lee” 

IF  Stepslow=l  THEN 

Kstart=(Kw(Stepslow+Extras)+Ks(Stepslow+Extras))/2 
Estart=(Ew(Stepslow+Extras)+Es(Stepslow+Extras))/2 
Knext=MIN(0.9999*Kl  ,0.99999*K2) 

ELSE 

IF  Stcpslow<Nmodes+l  THEN 

Kstart=(Kw(Stepslow+Extras)+Ks(Stepslow+Extras))/2 
Estart=^w(Stepslow+Extras)+Es(Stepslow+Extras))/2 
Kncxt=Ks(Stepslow-l ) 


ELSE 


Ksiart=l.E-4 

Estan=(Ew(Stepslow+Extras)+Es(Stepslow+Extras))A: 

Knext=Ks(Stepslow-l) 

END  IF 
END  IF 
END  SELECT 

PRINT 

Num_divs=INT(Onum_divs*(  1  +2*Stepsiow/(Nmodes+ 1 ))) 
Kstep=(Knext-Kstart)/Num_divs 

FOR  L=1  TO  Num_divs 
Num_iters=0 

Wip  1  ( 1 . 1  )=Kstan-^(L- 1  )*Kstep 

Wipl(2,l)=Estan 

IF  Wipl(l,l)oO  THEN 

REPEAT 
MAT  Wi=  Wipl 
GOSUB  Function 

IF  ABS(F(1,1))>100  OR  ABS(F(2,1))>100  THEN  Nextl 
GOSUB  Jacobian 
MAT  Wipl  =Jinv*F 

IF  Diffcase=l  THEN 

MAT  Wip  1=  Wip  1/(1 .0) 

ELSE 

MAT  Wipl=Wipl/(1.2) 

END  IF 

MAT  Wipl=  Wi-Wipl 
Num_iiers=Num_iters+ 1 

FOR  1=1  TO  Nmoden 

Tesi2=AB  S(Kn(I)-Wip  1(1,1)) 

IF  Tesi2<=0.0001  THEN  Next! 

NEXT  I 

IF  ABS(Wipl(l.l))>LE+5  THEN  Nextl 
IF  ABS(Wipl(2,l))>LE+5  THEN  Nextl 
F1=ABS(F(1,1)) 

F2=ABS(F(2.l)) 

UNTIL  Fl<=l.E-4  AND  F2<l.E-4  OR  Num_iters>=Max_iters 

IF  Num_iters>=Max_iters  THEN 
PRINT  “Max  Iterations” 

ELSE 

IF  Wip  1  ( 1 , 1  )<0  THEN  Nextl 
IF  Phil<0  Then  Nextl 
Nmodcn-*-Nmodcn+ 1 
MsNmoden 
Kn(M)=Wipl(l.l) 
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En(M)=Wip  1(2,1) 

SELECT  Type$ 

CASE  “111” 

PRINT  USING  Frmt2:”Kn(“.M,”)  =  “,Kn(M),”En(“,M,”)  =  “,En(M), 

R 1211,  R2311,Phil*  180/PI 
CASE  “lee” 

PRINT  USING  Frmtl;’-Knl(“,M,”)  =  “,Kn(M),”Enl(“,M,”)  =  “, 
En(M),Phil*  180/PI 
END  SELECT 
IF  Nmoden=l  THEN  Nextl 
GOTO  Nextss 
END  IF 
END  IF 
Nextl;! 

NEXTL 

Nextss:! 

NEXT  Stepslow 
GOTO  Sub_ex_end 

Function:! 

Km=Wi(l,l) 

Em=Wi(2,l) 

IF  Km=0  THEN  Nextl 

Kappa  1  =SQR(K  1  ''2-Km''2) 

Kappa2=SQR(K2''2-Km''2) 

Kappa2S=SQR(K2S'^2-Km''2) 

Signl=SGN(Kl'^2-Km''2) 

Sigr.2=SGN(K2'^2-Km''2) 

Sign2s=SGN(K2s''2-Km''2) 

IF  ABS(Km)>=ABS(K3)  THEN 
Kappa3=0 
ELSE 

Kappa3=SQR(K3''2-Kni'^2) 

END  IF 

Kappa  1  =Kappa  1  *Sign  1 
Kappa2=Kappa2  *  Sign2 
Kappa2s=Kappa2s*Sign2s 

Phil=ATN(Kappal/Ktn) 

Phi2=ATN(Kappa2/Kin) 

Rng=T/SIN(Phi2) 

D=Rhol/Rho2 
P=Kappa2/Kappa  1 
Z12=D*P 

Z23=(Rho2/Rho3)*(Kappa3/Kappa2) 

R1211=ABS((1-Z12)/(1+Z12)) 

R23U=ABS((  l-Z23)/(  1 +Z23)) 


Aa=2*  Kappa  1  *H+R  1  phase 
Bb=2*Kappa2*T 
Ab=Aa+Bb 
Saa=SIN(Aa) 

Sab=SIN(Ab) 

Caa=COS(Aa) 

Cab=COS(Ab) 

G=Km*Em 

G1=H/Kappal 

G2=T/Kappa2 

G12=G1+G2 

Gf=2*G*Gl 

Factl=^2*Km*Gl 

Fact2=-2*Km*G2 

Fact  1 2=Fact  1  +Fact2 

Ex2=-T*G/Kappa2 

Ex4=2*Ex2 

IF  ABS(Ex4)>10  OR  ABS(Gf)>10THEN  Nextl 

Y  l=Kin/Kappal 
>  2=Kin/Kappa2 
Y=1+Y1''2 
Y4=l+Y2''2 

B 1  =-2*  Alfa2*Rng*Freq/10000 
Fg=Gl*EXP(GO 

SELECT  Types 
CASE  “111” 

S=10^(B  1  )*(  1  -R  1211''2)*R2311*EXP(Ex4) 

S2=R  1 2U*R2311*  KV'IB  1  )*EXP(Ex4) 
Den=l+S2''2+2*S2*COS(Bb) 

F(  1 , 1  )=R  12ll*Saa+S*(Sab+Saa*S2)/Dcn 
F(2, 1  )=EXP(Gf>-R  1211*Caa^R  1  *S*(Cab+S2*Caa)/Dcn 
CASE  “lee” 

IF  ABS(Kin)>=ABS(K3s)  THEN 
Kappa3s=0 
ELSE 

Kappa3s=SQR(K3s''2-Km''2) 

END  IF 

Phi3=ATN(Kappa3/Km) 

Phi2s=ATN(Kappa2s/Km) 

Phi3s=ATN(Kappa3s^tn) 

Rngs=T/SIN(Phi2s) 

Cc=2*Kappa2s*T 

Dd=Kappa2*T+Kappa2s*T 

Ac=Aa+^ 

Sac=SIN(Ac) 

SadsSIN(Ad) 

Cac=COS(Ac) 

Cad=COS(Ad) 
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G2s=T/Kappa2s 

G12s=Gl+G2s 

Gs=2*GI+G2+G2s 

Ex2s=-T*G/Kappa2s 

Ex4s=2*Ex2s 

Ex=Ex2+Ex2s 

Y2s=Km/Kappa2s 

Y4s=l+Y2s'^2 

Cc  1  ^2*  Alfa2s*Rngs*Freq/10000 

D 1  ^(Alfa2*Rng+ AIfa2s*Rngs)*Freq/l 0000 

B12=l-(SIN(2*Phi2s))'^2*(l-Kappa2/Kappa2s) 

R121e=(-2l2+B12)/(zl2+B12) 

T12=2*D*COS(2*Phi2s)/(Z12+B  12) 

T2 1  =2*P*COS(2*Phi2s)/(Z  1 2+B 1 2) 

T1 2ps=-4*D*(COS(Phi2s)''2)*TAN(Phi2)/(Zl  2+B 1 2) 
T21sp^*P*(COS(Phi2s)'^2)*TAN(Phi2s)/(Z12+B12) 

R2  lden=SIN(Phi  1  )*((C2/C2s)*COS(2*Phi2s)^2+SIN(2*Phi2)*SIN(2*Phi2s)* 
(C2s/C2))+(C2*Cl*Rhol*SIN(Phi2))/(C2s*C2*Rho2) 
R21=((SIN(Phil)*((C2/C2s)*COS(2*Phi2s)''2-SIN(2*Phi2)*SIN(2*Phi2s)* 
(C2s/C2)HC2*C  1  *Rhol  *SIN(Phi2))/(C2s*C2*Rho2)))/R2  Iden 
R2 1  ss=((SIN(Phi  1  )*((C2/C2s)*COS(2*Phi2s)^2-SIN(2*Phi2)*SIN(2*Phi2s)* 
(C2s/C2))+(C2*Cl*Rhol*SIN(Phi2))/(C2s*C2*Rho2)))/R21den 

Tan2=Kappa2/Km 

Tan3=Kappa3/Km 

Tan2s=Kappa2s/Kin 

Tan3s=Kappa3s/Ktn 

Mu2s=Rho2*C2s''2 

Mu3s=Rho3*C3s''2 

Pl=(2*Mu2s+Mu3s*(Tan3s''2-l))*Tan2*Tan2s 

P2=(Mu2s*(Tan2s''2-l)-Mu3s*(Tan3s''2-l))*Tan2 

P3=(2*Mu3s+Mu2s*(Tan2s''2-l))*Tan3*Tan2s 

P4=2*(Mu3s-Mu2s)*Tan3*Tan2*Tan2s 

Q 1  =2*(Mu2s-Mu3s)*Tan3*Tan2*Tan2s 
Q2=(2*Mu3s+Mu2s*(Tan2s''2-l))*Tan3s*Tan2 
Q3=((-Mu2s*(Tan2s''2-l  ))+(Mu3s*(Tan3s^2-l  )))*Tan2s 
Q4=P1 


Denomin=((Pl+P3)*(Q2+Q4)-<P2+P4)*(Ql+Q3)) 
R23pp=(((Pl-P3)*(Q2+Q4)-(P2->-P4)*(Ql-Q3))/Denomin) 
R23ss=(((P2-P4)*(Ql  +Q3HP  l+P3)*(Q2-Q4))/Denomin) 
R23sp=(2*(P3*Ql-Pl  *Q3)/Dcnomin) 
R23ps=(2*(P4*Q2-P2*Q4)/Dcnomin) 

AsABS(Rl*R121c) 

P«ABS(R1*R23*T12*T21) 

C*ABS(R1  *R23ss*T12ps*T21sp) 
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D=ABS(R  I  *R23sp*T  1 2ps*T2 1 ) 

E=ABS(R  1  ♦R23ps*T12*T2 1  sp) 

B2=10'(Bl)*EXP(Ex4) 

Cc2=10^(Cc  1  )*EXP(Ex4s) 

D2=10^(Dl)*EXP(Ex) 

Pp=R2 1  *R23*  1 0^(B  1  )*EXP(Ex4) 

Ss=R2 1  ss*R23ss*  1 0'^CCc  1  )*EXP{Ex4s) 

Deno  1 = 1 -2*Pp*COS(Bb)-Pp''2 
Deno2=l-2*Ss*COS(Cc)-Ss^2 

FI  wat=A*EXP(-Gf)*Saa 

FI  sed=(B*  iO^(B  1  )/Deno  1  )*(EXP(-2*G*G  12)*Sab-Pp*EXP(-GO*Saa) 
Flscdsp=(D*10^(Dl)/Denol)*(EXP(-G*Gs)*Sad-Pp*EXP(-G*(2*Gl+G2s- 
G2))*SIN(Aa+Cc/2-Bb/2)) 

Flsedss=(C*10^(Ccl)/Deno2)*(EXP(-2*G*G12s)*Sac-Ss*EXP(-Gf)*SIN(Aa)) 
FI  sedps=(E*  lO'^CE  1  )/Deno2)*(EXP(-G*Gs)*Sad-Ss*EXP(-G*(2*G  1- 
G2s4G2))*SIN(Aa-Cc/2+Bb/2)) 

F2wat=A*Caa 

F2sed=(B*10^(Bl)/Denol)*(EXP(-2*G*G12)*Cab-Pp*EXP(-Gf)*Caa) 

FI  scdsp=(D*  10^(D1  )/Denol  )*(EXP(-G*Gs)*Cad-Pp*EXP(-<3*(2*Gl+G2s- 
G2))*COS(Aa+Cc/2-Bb/2)) 

FI  sedss=(C*  10^(Cc  1  )/Deno2)*(EXP(-2*G*G  1 2s)*Cac-Ss*EXP(-GO* 
COS(Aa)) 

Flsedps=(E*l(K'(El)/Deno2)*(EXP(-G*Gs)*Cad-Ss*EXP(-G*(2*Gl-- 

G2s+G2))*COS(Aa-Cc/2+Bb/2)) 


F(  1 , 1  )=F  1  wat+F  1  sed+F  1  sedss+F  1  sedsp+F  1  sedps 
F(2 , 1  )= 1-F2  wat*  EXP(-Gf)-F  1  sed-F  1  sedss-F  1  sedsp-F  1  sedps 
END  SELECT 
RETURN 


Jacobian:! 

SELECT  Type$ 

CASE“m” 

J(l,l)=Factl*R12U*Caa+S*(EXP(EX4)*(Cab*(Factl+Fact2>-Sab*2*Em*Y4)) 

J(  1  ^)=S*EXP(Ex4)*Sab*Fact2 

J(2, 1  )=EXP(Gf)*2*Em*G  1  *  Y^tR  1  ♦R  1211*Saa*Facil  +R1  *S*EXP(Ex4)* 
(Cab*2*Em*G2*Y4+Sab*(Factl+Fact2)) 

J(2.2)^EXP(Gf)*Fact  1-R 1  *S*EXP(Ex4)*Cab*Fact2 
CASE  “lee” 

Fact2s=^2*Km*G2s 
Fact  1 2=Fact  1  +Fact2 
Fact  1 2s=:Fact  1  +Fact2s 

Dfact=EXP(-G*Gs)*(Factl+Fact2/2+Fact2s/2) 

J(1 , 1  )=-A*EXP(-Gf)*Caa*Fact  1 +B*  10^(B  1  )*(EXP(-2*G*G  12)*Cab*Fact  12- 
Pp*EXP(-Gf)*Caa*Fact  1  )/Dcnol+C*  10^(Cc  1  )*(EXP(-2*G*G12s)*Cac* 
Factl2-Ss*EXP(-Gf)*Caa*Factl)/Dcno2+D*  K^HD  1  )*(Dfact*Cad-Pp* 
EXP(-G*(2*Gl-K32s^2))*COS(Aa+Cc/2-Bb/2)*(Factl+Fact2s/2- 


75 

Fact2/2)))/Denol+E*10^(Dl)*(Dfact*Cad-Ss*EXP(-G*(2*Gl-G2s+G2)) 

*COS(Aa-Cc/2+Bb/2)*(Factl-Fact2s/2+Fact2/2)))/Deno2 

J(l,2)=A*EXP(-G0*Saa*Factl+B*l(y'(Bl)*(EXP(-2*G*G12)*Sab*Factl2- 

Pp*EXP(-Gf)*Saa*Factl  )/Denol +C*  10^(Cc  1  )*(EXP(-2*G*G  1 2s)*Sac* 
Factl2-5s*EXP(-G0*Saa*Factl)/Deno2+D*l(y'(Dl)*(Dfact*Sad-Pp* 
EXP(-G*(2*Gl+G2s-G2))*SIN(Aa+Cc/2-Bb/2)*(Factl+Fact2s/2- 
Fact2/2)))/Denol +E*  10^(D1  )*(Dfact*Sad-Ss*EXP(-G*(2*G  1-G2S+G2)) 
*SIN(Aa-Cc/2+Bb/2)*(Factl-Fact2s/2+Fact2/2)))/Deno2 
J(2.1)=J(1.2) 

J(2.2)=-J(l,l) 

END  SELECT 
MATJinv=  INV(J) 

RETURN 

Sub_cx_end:! 

SUBEND 


SUB  Nat_log_e(E(*),Ln_e(*),Nmode) 

INTEGER  I 
FOR  1=1  TO  Nmode 

Ln_e(I)=LOG(ABS(E(I))) 
NEXT  I 
SUBEND 


SUB  Ln_e_vs_k_plot(K(*),E(*),Ln_e(*), Nmode, Color, MarkS) 

CSIZE  2,4 
PEN  Color 
LORG5 

FOR  I=Nmode  TO  1  STEP  -1 
MOVE  K(I),-Ln_e(I) 

LABEL  Marks 
NEXT  I 

SUBEND 


SUB  Incohmt(Kw(*),Ew(*),Ks(*),Es(*),Kn(*),En(*),Nmodew,Nmodes,Nmoden, 
Indbex,Extras) 

Formatl:  IMAGE  4X,56A.4D.D,X,2A 
COM  /Wavenumber/  W,  Kl,  K2,  K2s,  K3s 
INTEGER  I,R,Z 


Z=100 

R=10000 


Insumw=0 

lnsums=0 

lnsumn=0 
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FOR  I=Extras+l  TO  Nmodew 
Kappa  1  =SQR(K  1  ^2-Kw''2) 

Incow=((COS(Kappal*Z)*SIN(Kappal*Z)*EXP(-ABS(Ew(I))*R))/Kw(I)^1.5r2 
Insumw=Insumw+incow 
NEXT  I 

FOR  1=1  TO  Nmodes 

Kappal  =SQR(K  1^2-Ks''2) 

Incos=((COS(Kappal*Z)*SIN(Kappal*Z)*EXP(-ABS(Es(I))*R))/Ks(I)'^1.5)''2 
Insums=Insums+Incos 
NEXT  I 

FOR  1=1  TO  Nmoden 

Kappal  =SQR(K  1  ''2-Kn''2) 

Incon=((COS(Kappal*Z)*SlN(Kappal*Z)*EXP(-ABS(En(l))*R))/Kn(l)''1.5r2 

lnsunin=lnsumn+lncon 

NEXTl 

lnsumall=lnsuinw+lnsuins 
lndball=10*LGT(lnsumall) 
lndbex=  1 0*LGT(lnsumn) 

PRINT 

PRINT  USING  Fonnatl  “Incoherent  (Water  +  Sediment)  -  Incoherent  Exacts  ”,Indball- 
Indbex,“dB” 

SUBEND 


SUB  Fouraxes(  Color,  INTEGER  Diffcase) 

COM  /Parameterl/  Cl,  C2,  Rhol,  Rho2,  H,  T,  Freq,  Alfa2,  Rl,  Rlphase 
COM  AVavenumber/  W,  Kl,  K2,  K2s,  K3s 

DEG 
LORG5 
PEN  Color 

Kmax=INT(Freq/2)*  I  .E-2 
SELECT  Diffcase 
CASE  1 

VIEWPORT  20,65,47,72 
CASE  2 

VIEWPORT  67,1 12,47,72 
CASE  3 

VIEWPORT  67,1 12,21,46 
CASE  4 


VIEWPORT  20,65,21,46 
END  SELECT 

WINDOW  0,Kmax,0,7 
AXES  0,2,0, 1,1 
MOVE  0.05,1 
PEN-1 

DRAWKmax,! 

PEN  Color 
CSIZE  3,.4 

FOR  J=1  TO  10 

Ang=COS(10*(J-l))*Kl 
MOVE  Ang,1.25 
DRAW  Ang,l 
IF  J=1  THEN 
PEN-1 

MOVE  Ang+Kmax/35,1 
DRAWKmax,! 

PEN  Color 
END  IF 

IFDiffcase=3  OR  Diffcase=4  THEN 
IF  J=1  THEN 

MOVE  Ang,.3 
LABEL  “0” 

END  IF 
IF  J=4  THEN 

MOVE  Ang,.3 
LABEL  “30” 

END  IF 
IF  J=7  THEN 

MOVE  Ang„3 
LABEL  “60” 

END  IF 
IF  J=10THEN 

MOVE  Ang,.3 
LABEL  “90” 

END  IF 
END  IF 
NEXTJ 
SUBEND 


SUB  Label4axes(Color) 

VIEWPORT  0,131,0,100 
CSIZE  2.7 

LORG  5 

DEG 

PEN  Color 
MOVE  31,72 


78 


LABEL  “FINE  FLUID’ 

MOVE  75,72 
LABEL  “MEDIUM” 

MOVE  26,46 
LABEL  “HARD” 

MOVE  84,46 

LABEL  “SEMICONSOLIDATED” 

CSIZE3 
MOVE  62,18 

LABEL  “GRAZING  ANGLE” 

LDIR90 

MOVE  9,49 

LABEL  “-ln(ABS(En))” 

LDIRO 
CSIZE  3,.4 
MOVE  18,72 
LABEL ‘7” 

MOVE  18,51 
LABEL  “1” 

MOVE  18,46 
LABEL ‘7” 

MOVE  18,25 
LABEL  “1” 


SUBEND 


SUB  Son(K(*),E(*),Nmode) 

INTEGER  I,  N 
FOR  N=1  TO  Nmode-1 
FOR  I=N+1  TO  Nmode 
IF  K(N)>K(I)  THEN 
Junk=K(N) 
K(N)=K(I) 
K(I)=Junk 
Dummy=E(N) 
E(N)=E(I) 
E(I)=Junk 
END  IF 
NEXT  I 
NEXTN 
SUBEND 


